diff --git a/python_bindings/src/halide/CMakeLists.txt b/python_bindings/src/halide/CMakeLists.txt
index d1d6aaad00bc..5b85d545965e 100644
--- a/python_bindings/src/halide/CMakeLists.txt
+++ b/python_bindings/src/halide/CMakeLists.txt
@@ -34,6 +34,7 @@ list(TRANSFORM native_sources PREPEND "halide_/")
 set(python_sources
     __init__.py
     _generator_helpers.py
+    imageio.py
     )
 
 # It is technically still possible for a user to override the LIBRARY_OUTPUT_DIRECTORY by setting
diff --git a/python_bindings/src/halide/halide_/PyBuffer.cpp b/python_bindings/src/halide/halide_/PyBuffer.cpp
index 67f642413bb2..5f90d999c383 100644
--- a/python_bindings/src/halide/halide_/PyBuffer.cpp
+++ b/python_bindings/src/halide/halide_/PyBuffer.cpp
@@ -264,8 +264,8 @@ py::object buffer_setitem_operator(Buffer<> &buf, const std::vector<int> &pos, c
 class PyBuffer : public Buffer<> {
     py::buffer_info info;
 
-    PyBuffer(py::buffer_info &&info, const std::string &name)
-        : Buffer<>(pybufferinfo_to_halidebuffer(info), name),
+    PyBuffer(py::buffer_info &&info, const std::string &name, bool reverse_axes)
+        : Buffer<>(pybufferinfo_to_halidebuffer(info, reverse_axes), name),
           info(std::move(info)) {
     }
 
@@ -278,8 +278,8 @@ class PyBuffer : public Buffer<> {
         : Buffer<>(b), info() {
     }
 
-    PyBuffer(const py::buffer &buffer, const std::string &name)
-        : PyBuffer(buffer.request(/*writable*/ true), name) {
+    PyBuffer(const py::buffer &buffer, const std::string &name, bool reverse_axes)
+        : PyBuffer(buffer.request(/*writable*/ true), name, reverse_axes) {
         // Default to setting host-dirty on any PyBuffer we create from an existing py::buffer;
         // this allows (e.g.) code like
         //
@@ -299,6 +299,30 @@ class PyBuffer : public Buffer<> {
     ~PyBuffer() override = default;
 };
 
+py::buffer_info to_buffer_info(Buffer<> &b, bool reverse_axes = true) {
+    if (b.data() == nullptr) {
+        throw py::value_error("Cannot convert a Buffer<> with null host ptr to a Python buffer.");
+    }
+
+    const int d = b.dimensions();
+    const int bytes = b.type().bytes();
+    std::vector<Py_ssize_t> shape(d), strides(d);
+    for (int i = 0; i < d; i++) {
+        const int dst_axis = reverse_axes ? (d - i - 1) : i;
+        shape[dst_axis] = (Py_ssize_t)b.raw_buffer()->dim[i].extent;
+        strides[dst_axis] = (Py_ssize_t)b.raw_buffer()->dim[i].stride * (Py_ssize_t)bytes;
+    }
+
+    return py::buffer_info(
+        b.data(),                             // Pointer to buffer
+        bytes,                                // Size of one scalar
+        type_to_format_descriptor(b.type()),  // Python struct-style format descriptor
+        d,                                    // Number of dimensions
+        shape,                                // Buffer dimensions
+        strides                               // Strides (in bytes) for each index
+    );
+}
+
 }  // namespace
 
 void define_buffer(py::module &m) {
@@ -317,31 +341,12 @@ void define_buffer(py::module &m) {
             // Note that this allows us to convert a Buffer<> to any buffer-like object in Python;
             // most notably, we can convert to an ndarray by calling numpy.array()
             .def_buffer([](Buffer<> &b) -> py::buffer_info {
-                if (b.data() == nullptr) {
-                    throw py::value_error("Cannot convert a Buffer<> with null host ptr to a Python buffer.");
-                }
-
-                const int d = b.dimensions();
-                const int bytes = b.type().bytes();
-                std::vector<Py_ssize_t> shape, strides;
-                for (int i = 0; i < d; i++) {
-                    shape.push_back((Py_ssize_t)b.raw_buffer()->dim[i].extent);
-                    strides.push_back((Py_ssize_t)(b.raw_buffer()->dim[i].stride * bytes));
-                }
-
-                return py::buffer_info(
-                    b.data(),                             // Pointer to buffer
-                    bytes,                                // Size of one scalar
-                    type_to_format_descriptor(b.type()),  // Python struct-style format descriptor
-                    d,                                    // Number of dimensions
-                    shape,                                // Buffer dimensions
-                    strides                               // Strides (in bytes) for each index
-                );
+                return to_buffer_info(b, /*reverse_axes*/ true);
             })
 
             // This allows us to use any buffer-like python entity to create a Buffer<>
             // (most notably, an ndarray)
-            .def(py::init_alias<py::buffer, const std::string &>(), py::arg("buffer"), py::arg("name") = "")
+            .def(py::init_alias<py::buffer, const std::string &, bool>(), py::arg("buffer"), py::arg("name") = "", py::arg("reverse_axes") = true)
             .def(py::init_alias<>())
             .def(py::init_alias<const Buffer<> &>())
             .def(py::init([](Type type, const std::vector<int> &sizes, const std::string &name) -> Buffer<> {
@@ -404,6 +409,14 @@ void define_buffer(py::module &m) {
 
             .def("copy", &Buffer<>::copy)
             .def("copy_from", &Buffer<>::copy_from<void, Buffer<>::AnyDims>)
+            .def("reverse_axes", [](Buffer<> &b) -> Buffer<> {
+                const int d = b.dimensions();
+                std::vector<int> order;
+                for (int i = 0; i < b.dimensions(); i++) {
+                    order.push_back(d - i - 1);
+                }
+                return b.transposed(order);
+            })
 
             .def("add_dimension", (void(Buffer<>::*)()) & Buffer<>::add_dimension)
 
diff --git a/python_bindings/src/halide/halide_/PyBuffer.h b/python_bindings/src/halide/halide_/PyBuffer.h
index e675f6dd5d60..804d4b869804 100644
--- a/python_bindings/src/halide/halide_/PyBuffer.h
+++ b/python_bindings/src/halide/halide_/PyBuffer.h
@@ -15,7 +15,7 @@ py::object buffer_getitem_operator(Buffer<> &buf, const std::vector<int> &pos);
 template<typename T = void,
          int Dims = AnyDims,
          int InClassDimStorage = (Dims == AnyDims ? 4 : std::max(Dims, 1))>
-Halide::Runtime::Buffer<T, Dims, InClassDimStorage> pybufferinfo_to_halidebuffer(const py::buffer_info &info) {
+Halide::Runtime::Buffer<T, Dims, InClassDimStorage> pybufferinfo_to_halidebuffer(const py::buffer_info &info, bool reverse_axes) {
     const Type t = format_descriptor_to_type(info.format);
     halide_dimension_t *dims = (halide_dimension_t *)alloca(info.ndim * sizeof(halide_dimension_t));
     _halide_user_assert(dims);
@@ -23,7 +23,11 @@ Halide::Runtime::Buffer<T, Dims, InClassDimStorage> pybufferinfo_to_halidebuffer
         if (INT_MAX < info.shape[i] || INT_MAX < (info.strides[i] / t.bytes())) {
             throw py::value_error("Out of range dimensions in buffer conversion.");
         }
-        dims[i] = {0, (int32_t)info.shape[i], (int32_t)(info.strides[i] / t.bytes())};
+        // Halide's default indexing convention is col-major (the most rapidly varying index comes first);
+        // Numpy's default is row-major (most rapidly varying comes last).
+        // We usually want to reverse the order so that most-varying comes first.
+        const int dst_axis = reverse_axes ? (info.ndim - i - 1) : i;
+        dims[dst_axis] = {0, (int32_t)info.shape[i], (int32_t)(info.strides[i] / t.bytes())};
     }
     return Halide::Runtime::Buffer<T, Dims, InClassDimStorage>(t, info.ptr, (int)info.ndim, dims);
 }
@@ -31,8 +35,8 @@ Halide::Runtime::Buffer<T, Dims, InClassDimStorage> pybufferinfo_to_halidebuffer
 template<typename T = void,
          int Dims = AnyDims,
          int InClassDimStorage = (Dims == AnyDims ? 4 : std::max(Dims, 1))>
-Halide::Runtime::Buffer<T, Dims, InClassDimStorage> pybuffer_to_halidebuffer(const py::buffer &pyb, bool writable) {
-    return pybufferinfo_to_halidebuffer(pyb.request(writable));
+Halide::Runtime::Buffer<T, Dims, InClassDimStorage> pybuffer_to_halidebuffer(const py::buffer &pyb, bool writable, bool reverse_axes) {
+    return pybufferinfo_to_halidebuffer(pyb.request(writable), reverse_axes);
 }
 
 }  // namespace PythonBindings
diff --git a/python_bindings/src/halide/halide_/PyCallable.cpp b/python_bindings/src/halide/halide_/PyCallable.cpp
index 0d0229cad40e..b3aacec9b13f 100644
--- a/python_bindings/src/halide/halide_/PyCallable.cpp
+++ b/python_bindings/src/halide/halide_/PyCallable.cpp
@@ -92,7 +92,10 @@ class PyCallable {
                     argv[slot] = b.raw_buffer();
                 } else {
                     const bool writable = c_arg.is_output();
-                    buffers.buffers[slot] = pybuffer_to_halidebuffer<void, AnyDims, MaxFastDimensions>(cast_to<py::buffer>(value), writable);
+                    const bool reverse_axes = true;
+                    buffers.buffers[slot] =
+                        pybuffer_to_halidebuffer<void, AnyDims, MaxFastDimensions>(
+                            cast_to<py::buffer>(value), writable, reverse_axes);
                     argv[slot] = buffers.buffers[slot].raw_buffer();
                 }
                 cci[slot] = Callable::make_buffer_qcci();
diff --git a/python_bindings/src/halide/imageio.py b/python_bindings/src/halide/imageio.py
new file mode 100644
index 000000000000..ac927af4acc0
--- /dev/null
+++ b/python_bindings/src/halide/imageio.py
@@ -0,0 +1,73 @@
+try:
+    import imageio.v2 as imageio
+except:
+    import imageio
+import numpy
+
+
+def is_interleaved(im):
+    """If the given ndarray is 3-dimensional and appears to have an interleaved
+       layout, return True. Otherwise, return False."""
+
+    # Assume that 'interleaved' will only apply to channels <= 4
+    return im.ndim == 3 and im.strides[2] == 1 and im.shape[2] in [1, 2, 3, 4]
+
+
+def _as_interleaved(im):
+    """If the given ndarray is 3-dimensional and appears to be planar layout,
+       return a view that is in interleaved form, leaving the input unchanged.
+       Otherwise, return the image ndarray unchanged.
+       Note that this call must be used with care, as the returnee may or may
+       not be a copy."""
+    if im.ndim == 3 and not is_interleaved(im):
+        return numpy.moveaxis(im, 0, 2)
+    else:
+        return im
+
+
+def _as_planar(im):
+    """If the given ndarray is 3-dimensional and appears to be interleaved
+       layout, return a view that is in planar form, leaving the input
+       unchanged. Otherwise, return the image ndarray unchanged.
+       Note that this call must be used with care, as the returnee may or may
+       not be a copy."""
+    if is_interleaved(im):
+        return numpy.moveaxis(im, 2, 0)
+    else:
+        return im
+
+
+def copy_to_interleaved(im):
+    """If the given ndarray is 3-dimensional and appears to be planar
+       layout, return a copy that is in interleaved form. Otherwise, return
+       an unchanged copy of the input. Note that this call will always return
+       a copy, leaving the input unchanged."""
+    return _as_interleaved(im).copy()
+
+
+def copy_to_planar(im):
+    """If the given ndarray is 3-dimensional and appears to be interleaved
+       layout, return a copy that is in planar form. Otherwise, return
+       an unchanged copy of the input. Note that this call will always return
+       a copy, leaving the input unchanged."""
+    return _as_planar(im).copy()
+
+
+def imread(uri, format=None, **kwargs):
+    """halide.imageio.imread is a thin wrapper around imagio.imread,
+       except that for 3-dimensional images that appear to be interleaved,
+       the result is converted to a planar layout before returning."""
+    return copy_to_planar(imageio.imread(uri, format, **kwargs))
+
+
+def imwrite(uri, im, format=None, **kwargs):
+    """halide.imageio.imwrite is a thin wrapper around imagio.imwrite,
+       except that for 3-dimensional images that appear to be planar,
+       a temporary interleaved copy of the input is made, which is used for
+       writing."""
+
+    # We can use _as_interleaved() here to save a possible copy; since the
+    # caller will never see the possibly-a-copy value, there should be no
+    # risk of possibly-different behavior between cases that need converting
+    # and cases that don't.
+    imageio.imwrite(uri, _as_interleaved(im), format, **kwargs)
diff --git a/python_bindings/test/apps/bilateral_grid_shell.py b/python_bindings/test/apps/bilateral_grid_shell.py
index 87bb1d034250..880ea670e31d 100644
--- a/python_bindings/test/apps/bilateral_grid_shell.py
+++ b/python_bindings/test/apps/bilateral_grid_shell.py
@@ -6,7 +6,7 @@
 from bilateral_grid_Adams2019 import bilateral_grid_Adams2019
 from bilateral_grid_Li2018 import bilateral_grid_Li2018
 from bilateral_grid_Mullapudi2016 import bilateral_grid_Mullapudi2016
-import imageio
+import halide.imageio
 import numpy as np
 import os
 import sys
@@ -24,9 +24,7 @@ def main():
     r_sigma = float(sys.argv[3])
     timing_iterations = int(sys.argv[4])
 
-    assert os.path.exists(input_path), "Could not find %s" % input_path
-
-    input_buf_u8 = imageio.imread(input_path)
+    input_buf_u8 = halide.imageio.imread(input_path)
     assert input_buf_u8.dtype == np.uint8
     # Convert to float32
     input_buf = input_buf_u8.astype(np.float32)
@@ -54,7 +52,7 @@ def main():
 
     output_buf *= 255.0
     output_buf_u8 = output_buf.astype(np.uint8)
-    imageio.imsave(output_path, output_buf_u8)
+    halide.imageio.imwrite(output_path, output_buf_u8)
 
     print("Success!")
     sys.exit(0)
diff --git a/python_bindings/test/apps/blur.py b/python_bindings/test/apps/blur.py
index d89383166763..58700dc5f167 100644
--- a/python_bindings/test/apps/blur.py
+++ b/python_bindings/test/apps/blur.py
@@ -1,7 +1,7 @@
 import halide as hl
 
 import numpy as np
-import imageio
+import halide.imageio
 import os.path
 
 # Return the directory to look in for test images:
@@ -44,13 +44,10 @@ def get_blur(input):
 
 def get_input_data():
     image_path = os.path.join(apps_images_dir(), "rgb.png")
-    assert os.path.exists(image_path), \
-        "Could not find %s" % image_path
-    rgb_data = imageio.imread(image_path)
-    print("rgb_data", type(rgb_data), rgb_data.shape, rgb_data.dtype)
+    rgb_data = halide.imageio.imread(image_path)
 
-    grey_data = np.mean(rgb_data, axis=2, dtype=np.float32).astype(rgb_data.dtype)
-    input_data = np.copy(grey_data, order="F")
+    grey_data = np.mean(rgb_data, axis=0, dtype=np.float32).astype(rgb_data.dtype)
+    input_data = np.copy(grey_data)
 
     return input_data
 
@@ -65,20 +62,17 @@ def main():
     input_image = hl.Buffer(input_data)
     input.set(input_image)
 
-    output_data = np.empty(input_data.shape, dtype=input_data.dtype, order="F")
+    output_data = np.empty(input_data.shape, dtype=input_data.dtype)
     output_image = hl.Buffer(output_data)
 
-    print("input_image", input_image)
-    print("output_image", output_image)
-
     # do the actual computation
     blur.realize(output_image)
 
     # save results
     input_path = os.path.join(apps_output_dir(), "blur_input.png")
     output_path = os.path.join(apps_output_dir(), "blur_result.png")
-    imageio.imsave(input_path, input_data)
-    imageio.imsave(output_path, output_data)
+    halide.imageio.imwrite(input_path, input_data)
+    halide.imageio.imwrite(output_path, output_data)
     print("\nblur realized on output image.",
           "Result saved at", output_path,
           "( input data copy at", input_path, ")")
diff --git a/python_bindings/test/apps/erode.py b/python_bindings/test/apps/erode.py
index 95cbf9e6d39f..44f4815b5190 100644
--- a/python_bindings/test/apps/erode.py
+++ b/python_bindings/test/apps/erode.py
@@ -5,7 +5,7 @@
 import halide as hl
 
 import numpy as np
-import imageio
+import halide.imageio
 import os.path
 
 # Return the directory to look in for test images:
@@ -48,12 +48,10 @@ def get_erode(input):
 
 def get_input_data():
     image_path = os.path.join(apps_images_dir(), "rgb.png")
-    assert os.path.exists(image_path), \
-        "Could not find %s" % image_path
-    rgb_data = imageio.imread(image_path)
+    rgb_data = halide.imageio.imread(image_path)
     print("rgb_data", type(rgb_data), rgb_data.shape, rgb_data.dtype)
 
-    input_data = np.copy(rgb_data, order="F")
+    input_data = np.copy(rgb_data)
 
     return input_data
 
@@ -70,7 +68,7 @@ def main():
     input_image = hl.Buffer(input_data)
     input.set(input_image)
 
-    output_data = np.empty(input_data.shape, dtype=input_data.dtype, order="F")
+    output_data = np.empty(input_data.shape, dtype=input_data.dtype)
     output_image = hl.Buffer(output_data)
 
     print("input_image", input_image)
@@ -82,8 +80,8 @@ def main():
     # save results
     input_path = os.path.join(apps_output_dir(), "erode_input.png")
     output_path = os.path.join(apps_output_dir(), "erode_result.png")
-    imageio.imsave(input_path, input_data)
-    imageio.imsave(output_path, output_data)
+    halide.imageio.imwrite(input_path, input_data)
+    halide.imageio.imwrite(output_path, output_data)
     print("\nerode realized on output image.",
           "Result saved at", output_path,
           "( input data copy at", input_path, ")")
diff --git a/python_bindings/test/apps/interpolate.py b/python_bindings/test/apps/interpolate.py
index d6569d107368..5e23aa3020b4 100644
--- a/python_bindings/test/apps/interpolate.py
+++ b/python_bindings/test/apps/interpolate.py
@@ -5,7 +5,7 @@
 import halide as hl
 
 from datetime import datetime
-import imageio
+import halide.imageio
 import numpy as np
 import os.path
 
@@ -162,12 +162,10 @@ def get_interpolate(input, levels):
 
 def get_input_data():
     image_path = os.path.join(apps_images_dir(), "rgba.png")
-    assert os.path.exists(image_path), "Could not find %s" % image_path
-
-    rgba_data = imageio.imread(image_path)
+    rgba_data = halide.imageio.imread(image_path)
 
     # input data is in range [0, 1]
-    input_data = np.copy(rgba_data, order="F").astype(np.float32) / 255.0
+    input_data = np.copy(rgba_data).astype(np.float32) / 255.0
     return input_data
 
 
@@ -179,11 +177,11 @@ def main():
 
     # preparing input and output memory buffers (numpy ndarrays)
     input_data = get_input_data()
-    assert input_data.shape[2] == 4
     input_image = hl.Buffer(input_data)
+    assert input_image.channels() == 4
     input.set(input_image)
 
-    input_width, input_height = input_data.shape[:2]
+    input_width, input_height = input_image.width(), input_image.height()
 
     t0 = datetime.now()
     output_image = interpolate.realize([input_width, input_height, 3])
@@ -201,8 +199,8 @@ def main():
     # save results
     input_path = os.path.join(apps_output_dir(), "interpolate_input.png")
     output_path = os.path.join(apps_output_dir(), "interpolate_result.png")
-    imageio.imsave(input_path, input_data)
-    imageio.imsave(output_path, output_data)
+    halide.imageio.imwrite(input_path, input_data)
+    halide.imageio.imwrite(output_path, output_data)
 
     print()
     print('blur realized on output image. Result saved at {} (input data copy at {})'.format(output_path, input_path))
diff --git a/python_bindings/test/apps/local_laplacian.py b/python_bindings/test/apps/local_laplacian.py
index 99d9b4e38061..6a0a215f8059 100644
--- a/python_bindings/test/apps/local_laplacian.py
+++ b/python_bindings/test/apps/local_laplacian.py
@@ -5,7 +5,7 @@
 import halide as hl
 
 import numpy as np
-import imageio
+import halide.imageio
 import os.path
 
 # Return the directory to look in for test images:
@@ -188,12 +188,10 @@ def upsample2D(f):
 
 def get_input_data():
     image_path = os.path.join(apps_images_dir(), "rgb.png")
-    assert os.path.exists(image_path), "Could not find {}".format(image_path)
-
-    rgb_data = imageio.imread(image_path)
+    rgb_data = halide.imageio.imread(image_path)
 
     # input data is in range [0, 256*256]
-    input_data = rgb_data.astype(np.uint16, order="F") << 8
+    input_data = rgb_data.astype(np.uint16) << 8
     return input_data
 
 
@@ -208,7 +206,7 @@ def filter_test_image(local_laplacian, input):
     output_data = np.empty_like(input_data)
 
     # do the actual computation
-    input_width, input_height = input_data.shape[:2]
+    input_width, input_height = input_image.width(), input_image.height()
     output_image = local_laplacian.realize([input_width, input_height, 3])
     output_data = np.asanyarray(output_image)
 
@@ -220,8 +218,8 @@ def filter_test_image(local_laplacian, input):
     input_path = os.path.join(apps_output_dir(), "local_laplacian_input.png")
     output_path = os.path.join(apps_output_dir(), "local_laplacian.png")
 
-    imageio.imsave(input_path, input_data)
-    imageio.imsave(output_path, output_data)
+    halide.imageio.imwrite(input_path, input_data)
+    halide.imageio.imwrite(output_path, output_data)
 
     print()
     print("local_laplacian realized on output_image.")
diff --git a/python_bindings/test/correctness/buffer.py b/python_bindings/test/correctness/buffer.py
index 10196300eb6e..39ba655a18cd 100644
--- a/python_bindings/test/correctness/buffer.py
+++ b/python_bindings/test/correctness/buffer.py
@@ -4,83 +4,131 @@
 import gc
 import sys
 
-def test_ndarray_to_buffer():
+def test_ndarray_to_buffer(reverse_axes = True):
     a0 = np.ones((200, 300), dtype=np.int32)
 
     # Buffer always shares data (when possible) by default,
     # and maintains the shape of the data source. (note that
     # the ndarray is col-major by default!)
-    b0 = hl.Buffer(a0, "float32_test_buffer")
+    b0 = hl.Buffer(a0, "float32_test_buffer", reverse_axes)
     assert b0.type() == hl.Int(32)
     assert b0.name() == "float32_test_buffer"
     assert b0.all_equal(1)
 
-    assert b0.dim(0).min() == 0
-    assert b0.dim(0).max() == 199
-    assert b0.dim(0).extent() == 200
-    assert b0.dim(0).stride() == 300
-
-    assert b0.dim(1).min() == 0
-    assert b0.dim(1).max() == 299
-    assert b0.dim(1).extent() == 300
-    assert b0.dim(1).stride() == 1
-
-    a0[12, 34] = 56
-    assert b0[12, 34] == 56
-
-    b0[56, 34] = 12
-    assert a0[56, 34] == 12
-
-
-def test_buffer_to_ndarray():
-    buf = hl.Buffer(hl.Int(16), [4, 4])
-    assert buf.type() == hl.Int(16)
-    buf.fill(0)
-    buf[1, 2] = 42
-    assert buf[1, 2] == 42
+    if reverse_axes:
+        assert b0.dim(0).min() == 0
+        assert b0.dim(0).max() == 299
+        assert b0.dim(0).extent() == 300
+        assert b0.dim(0).stride() == 1
+
+        assert b0.dim(1).min() == 0
+        assert b0.dim(1).max() == 199
+        assert b0.dim(1).extent() == 200
+        assert b0.dim(1).stride() == 300
+
+        a0[12, 34] = 56
+        assert b0[34, 12] == 56
+
+        b0[56, 34] = 12
+        assert a0[34, 56] == 12
+    else:
+        assert b0.dim(0).min() == 0
+        assert b0.dim(0).max() == 199
+        assert b0.dim(0).extent() == 200
+        assert b0.dim(0).stride() == 300
+
+        assert b0.dim(1).min() == 0
+        assert b0.dim(1).max() == 299
+        assert b0.dim(1).extent() == 300
+        assert b0.dim(1).stride() == 1
+
+        a0[12, 34] = 56
+        assert b0[12, 34] == 56
+
+        b0[56, 34] = 12
+        assert a0[56, 34] == 12
+
+
+def test_buffer_to_ndarray(reverse_axes = True):
+    buf0 = hl.Buffer(hl.Int(16), [4, 6])
+    assert buf0.type() == hl.Int(16)
+    buf0.fill(0)
+    buf0[1, 2] = 42
+    assert buf0[1, 2] == 42
+
+    # This is subtle: the default behavior when converting
+    # a Buffer to an np.array (or ndarray, etc) is to reverse the
+    # order of the axes, since Halide prefers column-major and
+    # the rest of Python prefers row-major. By calling reverse_axes()
+    # before that conversion, we end up doing a *double* reverse, i.e,
+    # not reversing at all. So the 'not' here is correct.
+    buf = buf0.reverse_axes() if not reverse_axes else buf0
 
     # Should share storage with buf
     array_shared = np.array(buf, copy = False)
-    assert array_shared.shape == (4, 4)
     assert array_shared.dtype == np.int16
-    assert array_shared[1, 2] == 42
+    if reverse_axes:
+        assert array_shared.shape == (6, 4)
+        assert array_shared[2, 1] == 42
+    else:
+        assert array_shared.shape == (4, 6)
+        assert array_shared[1, 2] == 42
 
     # Should *not* share storage with buf
     array_copied = np.array(buf, copy = True)
-    assert array_copied.shape == (4, 4)
     assert array_copied.dtype == np.int16
-    assert array_copied[1, 2] == 42
-
-    buf[1, 2] = 3
-    assert array_shared[1, 2] == 3
-    assert array_copied[1, 2] == 42
+    if reverse_axes:
+        assert array_copied.shape == (6, 4)
+        assert array_copied[2, 1] == 42
+    else:
+        assert array_copied.shape == (4, 6)
+        assert array_copied[1, 2] == 42
+
+    # Should affect array_shared but not array_copied
+    buf0[1, 2] = 3
+    if reverse_axes:
+        assert array_shared[2, 1] == 3
+        assert array_copied[2, 1] == 42
+    else:
+        assert array_shared[1, 2] == 3
+        assert array_copied[1, 2] == 42
 
     # Ensure that Buffers that have nonzero mins get converted correctly,
     # since the Python Buffer Protocol doesn't have the 'min' concept
-    cropped = buf.copy()
-    cropped.crop(dimension = 0, min = 1, extent = 2)
+    cropped_buf0 = buf0.copy()
+    cropped_buf0.crop(dimension = 0, min = 1, extent = 2)
+    cropped_buf = cropped_buf0.reverse_axes() if not reverse_axes else cropped_buf0
 
     # Should share storage with cropped (and buf)
-    cropped_array_shared = np.array(cropped, copy = False)
-    assert cropped_array_shared.shape == (2, 4)
+    cropped_array_shared = np.array(cropped_buf, copy = False)
     assert cropped_array_shared.dtype == np.int16
-    assert cropped_array_shared[0, 2] == 3
+    if reverse_axes:
+        assert cropped_array_shared.shape == (6, 2)
+        assert cropped_array_shared[2, 0] == 3
+    else:
+        assert cropped_array_shared.shape == (2, 6)
+        assert cropped_array_shared[0, 2] == 3
 
     # Should *not* share storage with anything
-    cropped_array_copied = np.array(cropped, copy = True)
-    assert cropped_array_copied.shape == (2, 4)
+    cropped_array_copied = np.array(cropped_buf, copy = True)
     assert cropped_array_copied.dtype == np.int16
-    assert cropped_array_copied[0, 2] == 3
-
-    cropped[1, 2] = 5
-
-    assert buf[1, 2] == 3
-    assert array_shared[1, 2] == 3
-    assert array_copied[1, 2] == 42
-
-    assert cropped[1, 2] == 5
-    assert cropped_array_shared[0, 2] == 5
-    assert cropped_array_copied[0, 2] == 3
+    if reverse_axes:
+        assert cropped_array_copied.shape == (6, 2)
+        assert cropped_array_copied[2, 0] == 3
+    else:
+        assert cropped_array_copied.shape == (2, 6)
+        assert cropped_array_copied[0, 2] == 3
+
+    cropped_buf0[1, 2] = 5
+    assert cropped_buf0[1, 2] == 5
+    if reverse_axes:
+        assert cropped_buf[1, 2] == 5
+        assert cropped_array_shared[2, 0] == 5
+        assert cropped_array_copied[2, 0] == 3
+    else:
+        assert cropped_buf[2, 1] == 5
+        assert cropped_array_shared[0, 2] == 5
+        assert cropped_array_copied[0, 2] == 3
 
 
 def _assert_fn(e):
@@ -169,8 +217,9 @@ def test_make_interleaved():
     assert b.dim(2).stride() == 1
 
     a = np.array(b, copy = False)
-    assert a.shape == (w, h, c)
-    assert a.strides == (c, w*c, 1)
+    # NumPy shape order is opposite that of Halide shape order
+    assert a.shape == (c, h, w)
+    assert a.strides == (1, w*c, c)
     assert a.dtype == np.uint8
 
 def test_interleaved_ndarray():
@@ -188,16 +237,16 @@ def test_interleaved_ndarray():
     assert b.type() == hl.UInt(8)
 
     assert b.dim(0).min() == 0
-    assert b.dim(0).extent() == w
-    assert b.dim(0).stride() == c
+    assert b.dim(0).extent() == c
+    assert b.dim(0).stride() == 1
 
     assert b.dim(1).min() == 0
     assert b.dim(1).extent() == h
     assert b.dim(1).stride() == w * c
 
     assert b.dim(2).min() == 0
-    assert b.dim(2).extent() == c
-    assert b.dim(2).stride() == 1
+    assert b.dim(2).extent() == w
+    assert b.dim(2).stride() == c
 
 def test_reorder():
     W = 7
@@ -288,8 +337,10 @@ def test_scalar_buffers():
 if __name__ == "__main__":
     test_make_interleaved()
     test_interleaved_ndarray()
-    test_ndarray_to_buffer()
-    test_buffer_to_ndarray()
+    test_ndarray_to_buffer(reverse_axes = True)
+    test_ndarray_to_buffer(reverse_axes = False)
+    test_buffer_to_ndarray(reverse_axes = True)
+    test_buffer_to_ndarray(reverse_axes = False)
     test_for_each_element()
     test_fill_all_equal()
     test_bufferinfo_sharing()
diff --git a/python_bindings/tutorial/lesson_02_input_image.py b/python_bindings/tutorial/lesson_02_input_image.py
index e24ca547237b..f132ff75ba1b 100644
--- a/python_bindings/tutorial/lesson_02_input_image.py
+++ b/python_bindings/tutorial/lesson_02_input_image.py
@@ -8,7 +8,7 @@
 
 import halide as hl
 import numpy as np
-import imageio
+import halide.imageio
 import os.path
 
 
@@ -21,7 +21,7 @@ def main():
     image_path = os.path.join(os.path.dirname(__file__), "../../tutorial/images/rgb.png")
 
     # We create a hl.Buffer object to wrap the numpy array
-    input = hl.Buffer(imageio.imread(image_path))
+    input = hl.Buffer(halide.imageio.imread(image_path))
     assert input.type() == hl.UInt(8)
 
     # Next we define our hl.Func object that represents our one pipeline
@@ -87,7 +87,7 @@ def main():
 
     # Save the output for inspection. It should look like a bright parrot.
     # python3-imageio versions <2.5 expect a numpy array
-    imageio.imsave("brighter.png", np.asanyarray(output_image))
+    halide.imageio.imwrite("brighter.png", np.asanyarray(output_image))
     print("Created brighter.png result file.")
 
     print("Success!")
diff --git a/python_bindings/tutorial/lesson_07_multi_stage_pipelines.py b/python_bindings/tutorial/lesson_07_multi_stage_pipelines.py
index db04954f0856..1800d0ebbaf9 100644
--- a/python_bindings/tutorial/lesson_07_multi_stage_pipelines.py
+++ b/python_bindings/tutorial/lesson_07_multi_stage_pipelines.py
@@ -10,7 +10,7 @@
 
 import halide as hl
 
-import imageio
+import halide.imageio
 import numpy as np
 import os.path
 
@@ -25,7 +25,7 @@ def main():
     # first horizontally, and then vertically.
     if True:
         # Take a color 8-bit input
-        input = hl.Buffer(imageio.imread(image_path))
+        input = hl.Buffer(halide.imageio.imread(image_path))
         assert input.type() == hl.UInt(8)
 
         # Upgrade it to 16-bit, so we can do math without it overflowing.
@@ -84,7 +84,7 @@ def main():
         # shorter than the input image.
 
         # python3-imageio versions <2.5 expect a numpy array
-        imageio.imsave("blurry_parrot_1.png", np.asanyarray(result))
+        halide.imageio.imwrite("blurry_parrot_1.png", np.asanyarray(result))
         print("Created blurry_parrot_1.png")
 
         # This is usually the fastest way to deal with boundaries:
@@ -94,7 +94,7 @@ def main():
     # The same pipeline, with a boundary condition on the input.
     if True:
         # Take a color 8-bit input
-        input = hl.Buffer(imageio.imread(image_path))
+        input = hl.Buffer(halide.imageio.imread(image_path))
         assert input.type() == hl.UInt(8)
 
         # This time, we'll wrap the input in a hl.Func that prevents
@@ -155,7 +155,7 @@ def main():
         # input.
 
         # python3-imageio versions <2.5 expect a numpy array
-        imageio.imsave("blurry_parrot_2.png", np.asanyarray(result))
+        halide.imageio.imwrite("blurry_parrot_2.png", np.asanyarray(result))
         print("Created blurry_parrot_2.png")
 
     print("Success!")
diff --git a/python_bindings/tutorial/lesson_09_update_definitions.py b/python_bindings/tutorial/lesson_09_update_definitions.py
index db93b7674f05..545f4424ddab 100644
--- a/python_bindings/tutorial/lesson_09_update_definitions.py
+++ b/python_bindings/tutorial/lesson_09_update_definitions.py
@@ -12,7 +12,7 @@
 
 import halide as hl
 
-import imageio
+import halide.imageio
 import numpy as np
 import os.path
 
@@ -23,10 +23,10 @@ def main():
 
     # Load a grayscale image to use as an input.
     image_path = os.path.join(os.path.dirname(__file__), "../../tutorial/images/gray.png")
-    input_data = imageio.imread(image_path)
+    input_data = halide.imageio.imread(image_path)
     if True:
          # making the image smaller to go faster
-        input_data = input_data[:160, :150]
+        input_data = input_data[:150, :160]
     assert input_data.dtype == np.uint8
     input = hl.Buffer(input_data)
 
@@ -193,7 +193,7 @@ def main():
 
         for r_y in range(input.height()):
             for r_x in range(input.width()):
-                py_result[input_data[r_x, r_y]] += 1
+                py_result[input_data[r_y, r_x]] += 1
 
         # Check the answers agree:
         for xx in range(256):
diff --git a/python_bindings/tutorial/lesson_12_using_the_gpu.py b/python_bindings/tutorial/lesson_12_using_the_gpu.py
index 913c55e2cc9e..e0c22615bb93 100644
--- a/python_bindings/tutorial/lesson_12_using_the_gpu.py
+++ b/python_bindings/tutorial/lesson_12_using_the_gpu.py
@@ -10,7 +10,7 @@
 
 import halide as hl
 
-import imageio
+import halide.imageio
 import os.path
 import struct
 
@@ -240,7 +240,7 @@ def test_correctness(self, reference_output):
 def main():
     # Load an input image.
     image_path = os.path.join(os.path.dirname(__file__), "../../tutorial/images/rgb.png")
-    input = hl.Buffer(imageio.imread(image_path))
+    input = hl.Buffer(halide.imageio.imread(image_path))
 
     # Allocated an image that will store the correct output
     reference_output = hl.Buffer(hl.UInt(8), [input.width(), input.height(), input.channels()])