python_cpp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "![Translation](trans.png)\n",
    "\n",
    "# Calling C++ code from Python\n",
    "\n",
    "## Problem\n",
    "\n",
    " - We have some existing C++ code which operates on array/image data\n",
    " - We want to call it from Python\n",
    " - We want to use Numpy arrays to pass input and receive output\n",
    " - **Ideally, want to avoid too much copying of large data**\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution I will present\n",
    "\n",
    " - Build a **Cython** extension\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alternative solutions I will mention briefly\n",
    "\n",
    " - Create a pure-C API and use `ctypes`\n",
    " - Wrapper-generators (e.g. `swig`)\n",
    " - Wrapping a command line tool\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sample C++ code\n",
    "    #include \"newimage/newimageall.h\"\n",
    "\n",
    "    void process_volume(NEWIMAGE::volume4D<float> &invol)\n",
    "    {\n",
    "        // Do some clever stuff\n",
    "        invol.binarise(0.5);\n",
    "    }\n",
    "\n",
    "    int main(int argc, char **argv)\n",
    "    {\n",
    "        char *input_file = argv[1];\n",
    "        char *output_file = argv[2];\n",
    "\n",
    "        NEWIMAGE::volume4D<float> invol;\n",
    "        read_volume4D(invol, input_file);\n",
    "\n",
    "        process_volume(invol);\n",
    "        save_volume4D(invol, output_file);\n",
    "    }\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    \n",
    "## First provide an entry point using C++ native types\n",
    "\n",
    "    #include <vector>\n",
    "    #include <iostream>\n",
    "    \n",
    "    std::vector<float> process_vectors(std::vector<float> &input, int nx, int ny, int nz, int nt)\n",
    "    {\n",
    "        // This is just so we can see if the data has been copied\n",
    "        std::cerr << \"In C++ the input vector starts at address \" << &input[0] << std::endl;\n",
    "        \n",
    "        // Here we ought to check that nx, ny, nz, nt is consistent with overall length of input\n",
    "        \n",
    "        // Create a volume4D using an existing data buffer\n",
    "        // when we do this, NEWIMAGE will not try to delete the data buffer\n",
    "        NEWIMAGE::volume4D<float> invol(nx, ny, nz, nt, &input[0]);\n",
    "\n",
    "        // Do our processing step\n",
    "        process_volume(invol);\n",
    "        \n",
    "        // Input data has been modified, so return it directly\n",
    "        return input;\n",
    "    }\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Array ordering\n",
    "\n",
    "![Ordering](abc.png)\n",
    "\n",
    "If `input` is a 4D image, it's pretty clear that the first element is the voxel with co-ordinates `(0, 0, 0, 0)`\n",
    "\n",
    "But what is the next element?\n",
    "\n",
    "Is it voxel `(1, 0, 0, 0)`?\n",
    "\n",
    "Or `(0, 0, 0, 1)`?\n",
    "\n",
    "If the *first* axis is the one which varies fastest, we are using **Column-Major** ordering\n",
    "If the *last* axis is the one which varies fastest, we are using **Row-Major** ordering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## So, which is the standard?\n",
    "\n",
    "| Row-major            | Column-major | \n",
    "| ------------         | ---------    | \n",
    "| C/C++ native arrays  | Fortran      | \n",
    "| Python/Numpy default | Matlab       |\n",
    "| SAS                  | FSL NEWIMAGE | \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Here, we will need to make sure our Numpy arrays are passed as 1-dimensional float arrays in Column-major order to match NEWIMAGE\n",
    "\n",
    "    data.flatten(order='F').astype(np.float32)\n",
    "    \n",
    " - `'F'` stands for 'Fortran order'\n",
    " - A C++ `float` is *almost* guaranteed to be 32 bits\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "![Cython](cython.png)\n",
    "\n",
    "# Cython extension\n",
    "\n",
    "## First, the Cython wrapper\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# my_analysis_wrapper.pyx\n",
    "\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "\n",
    "from libcpp.vector cimport vector\n",
    "\n",
    "cdef extern from \"my_analysis.h\":\n",
    "    vector[float] process_vector(vector[float] &, int, int, int, int)\n",
    "    \n",
    "def process_using_vectors(data):\n",
    "    # Save the dimensions of the data because we're going to flatten it to 1D array\n",
    "    # Should be checking the dimensions at this point!\n",
    "    nx, ny, nz, nt = data.shape\n",
    "\n",
    "    # Convert data to 1D in Column-major (Fortran) order\n",
    "    # This always copies the data\n",
    "    data = data.flatten(order='F').astype(np.float32)\n",
    "\n",
    "    # This line is just so we can see if the data is being copied\n",
    "    print(\"In python the input data starts at %X\" % data.__array_interface__['data'][0])\n",
    "\n",
    "    # Call the C++ code\n",
    "    output = process_vectors(data, nx, ny, nz, nt)\n",
    "\n",
    "    # Output is a 1D array in Fortran order - turn it back into a multidimensional array\n",
    "    # This should not copy the data\n",
    "    output = np.reshape(output, [nx, ny, nz, nt], order='F')\n",
    "    print(\"In python the reshaped data starts at %X\" % output.__array_interface__['data'][0])\n",
    "    \n",
    "    return output\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        \n",
    "## Next, build the extension\n",
    "\n",
    "This would normally go in `setup.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import numpy\n",
    "\n",
    "from setuptools import setup\n",
    "from Cython.Build import cythonize\n",
    "from setuptools.extension import Extension\n",
    "\n",
    "# My Cython extension\n",
    "fsldir = os.environ[\"FSLDIR\"]\n",
    "\n",
    "ext = Extension(\"my_analysis_wrapper\",\n",
    "                 sources=['my_analysis_wrapper.pyx',\n",
    "                          'my_analysis.cpp'],\n",
    "                 language=\"c++\",\n",
    "                 include_dirs=[\".\", numpy.get_include(), \n",
    "                               os.path.join(fsldir, \"include\"), \n",
    "                               os.path.join(fsldir, \"extras/include\"), \n",
    "                               os.path.join(fsldir, \"extras/include/newmat\")], \n",
    "                 libraries=['newimage', 'miscmaths', 'fslio', 'niftiio', 'newmat', 'znz', \"zlib\"],\n",
    "                 library_dirs=[os.path.join(fsldir, \"lib\"), os.path.join(fsldir, \"extras/lib\")])\n",
    "\n",
    "# setup parameters\n",
    "setup(name='my_app',\n",
    "      description='My Python application which calls C++',\n",
    "      ext_modules=cythonize(ext))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running build_ext\n",
      "copying build\\lib.win-amd64-2.7\\my_analysis_wrapper.pyd -> \n"
     ]
    }
   ],
   "source": [
    "%run setup.py build_ext"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD8CAYAAABaQGkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAADKRJREFUeJzt3W1sXvV5x/HfL3Yc7DQkTbulysPq\nwChrSFeBPMaDxgR0E5SMrGonkYpqY9WiTQVS2q2ibBKa9mJvECqqKjqPPkhrVMbSjFUdg9LSFqFp\nVpyEKhiXNktYHoGUpCTQ0Djk2gt7UsqI72P8/3PsS9+PhBSbw8UlJ9+c2/d9fG5HhADkNKftBQDU\nQ+BAYgQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJNZdY2jPnN7o7VpQfO6C839RfKYkdelU8Zkv7l9U\nfKYkvTavylh1v1rnisbXelxl7pIlR4rPPPDKwuIzJanr5fLn0RPHDuvk8Vc6fnGrBN7btUCXLv5I\n8blXPrCr+ExJWjDn1eIzN/71muIzJemllV1V5i5+ZqzK3JfePbfK3L/c8M/FZ9659Q+Kz5SkBU/0\nFZ/5kwfubnQcD9GBxAgcSIzAgcQIHEiMwIHECBxIrFHgtq+x/YztnbZvr70UgDI6Bm67S9IXJF0r\naZWkdbZX1V4MwPQ1OYNfLGlnROyKiBOS7pe0tu5aAEpoEvgySXtP+3jfxOd+ie31todtD584dbzU\nfgCmoUngb3S96/+7cDkiBiNiICIGeub0Tn8zANPWJPB9klac9vFySQfqrAOgpCaBb5F0nu2Vtnsk\n3SDpm3XXAlBCx58mi4iTtm+W9IikLklfjoiR6psBmLZGPy4aEQ9JeqjyLgAK40o2IDECBxIjcCAx\nAgcSI3AgsSo3XdTcbsWSdxQf+/TL5W+OKEkjX1xdfObLq+vcTfS3r9tRZe6WB99XZe7Y/Dp3a733\njvI39Tx/68HiMyVpbGn5r8Hunze7EzBncCAxAgcSI3AgMQIHEiNwIDECBxIjcCAxAgcSI3AgMQIH\nEiNwIDECBxIjcCAxAgcSI3AgMQIHEiNwIDECBxIjcCAxAgcSI3AgsTp3VZWk7vJ/d7zwh33FZ0rS\ndY/8oPjMB5/9zeIzJenx/7qgytwf3nx3lbm/9Y+fqjJ3/v7jxWf++3/WedPccx+7qfjMX/xNs+M4\ngwOJETiQGIEDiRE4kBiBA4kROJBYx8Btr7D9Pdujtkdsb3grFgMwfU1eBz8p6dMRsc32AklbbT8a\nEU9X3g3ANHU8g0fEwYjYNvHrY5JGJS2rvRiA6ZvS9+C2+yVdKGmoxjIAymocuO23SfqGpE9GxNE3\n+PfrbQ/bHj5x8ucldwTwJjUK3PZcjce9MSI2v9ExETEYEQMRMdDTXeeacQBT0+RZdEv6kqTRiKjz\nEwkAqmhyBr9c0sckXWX7yYl/Plh5LwAFdHyZLCKekOS3YBcAhXElG5AYgQOJETiQGIEDiRE4kFiV\nmy6+urhLO9ctLD73u+sGi8+UpFt3f7j4zL94z+PFZ0rS379Y5xXK9z36iSpz+4fGqsz98Z/OKz7z\nypG1xWdK0rJ/mVt85qEjzV7Y4gwOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGIEDiRG4EBiBA4k\nRuBAYgQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGIEDiRW5a6qPfPHtPyiA8Xn3nPoiuIzJWnX\nv51bfOaDG6P4TEla+rtVfst01vryv1+StHvdO6vMXflP5d8u76wf1nlf+1NHthWfOedks105gwOJ\nETiQGIEDiRE4kBiBA4kROJAYgQOJNQ7cdpft7ba/VXMhAOVM5Qy+QdJorUUAlNcocNvLJV0n6b66\n6wAoqekZ/HOSPiPp1JkOsL3e9rDt4bGf1bnkD8DUdAzc9hpJL0TE1smOi4jBiBiIiIG5i/qKLQjg\nzWtyBr9c0vW2n5V0v6SrbH+t6lYAiugYeER8NiKWR0S/pBskPRYRN1bfDMC08To4kNiUfrg4Ir4v\n6ftVNgFQHGdwIDECBxIjcCAxAgcSI3AgsSq36Bw72aUDhxcWn3vBoueKz5SkFf+6v/jMf9iyufhM\nSfrgtj+rMvfs639aZe7cO5ZVmfvcnx8rPrP30V8vPlOSeo6Vv8Pua//xRKPjOIMDiRE4kBiBA4kR\nOJAYgQOJETiQGIEDiRE4kBiBA4kROJAYgQOJETiQGIEDiRE4kBiBA4kROJAYgQOJETiQGIEDiRE4\nkBiBA4k5ovwdHxfOWxKXveujxec+//sris+UpMO/c6L4zPO+OFZ8piS9uHp+lbmrPj5SZe7zlx6t\nMvfjP95dfObQsXOLz5Sk31tY/mu7Ye1/6yc7jrvTcZzBgcQIHEiMwIHECBxIjMCBxAgcSKxR4LYX\n2d5k+0e2R21fWnsxANPX9N1F75H0cER8xHaPpL6KOwEopGPgts+WdIWkP5GkiDghqfyVIQCKa/IQ\n/RxJhyR9xfZ22/fZrnM5FYCimgTeLekiSfdGxIWSXpF0++sPsr3e9rDt4ROvHS+8JoA3o0ng+yTt\ni4ihiY83aTz4XxIRgxExEBEDPV29JXcE8CZ1DDwinpO01/b5E5+6WtLTVbcCUETTZ9FvkbRx4hn0\nXZJuqrcSgFIaBR4RT0oaqLwLgMK4kg1IjMCBxAgcSIzAgcQIHEiMwIHEmr4OPjX9UgyeKj72Vz9U\n5/qaJd85u/jMXXctKj5Tkvr/7kiVuUP9F1SZu7J7S5W5f7tjTfGZffPq3An34WffW3zmvuODjY7j\nDA4kRuBAYgQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGIE\nDiRG4EBiBA4kRuBAYlVuuvhr8w7r8+c8UHxu/0hf8ZmS9NWjS4vPvH/9tcVnStLO27uqzP3oBY9X\nmfvwM1dUmfvq3ig+c9sffb74TEn60PuvKT7zwJFmN4jkDA4kRuBAYgQOJEbgQGIEDiRG4EBiBA4k\n1ihw27fZHrH9lO2v2z6r9mIApq9j4LaXSbpV0kBErJbUJemG2osBmL6mD9G7JfXa7pbUJ+lAvZUA\nlNIx8IjYL+kuSXskHZT0UkR8+/XH2V5ve9j28OHD5d8bHMDUNXmI/nZJayWtlLRU0nzbN77+uIgY\njIiBiBhYvJjn7oCZoEmJH5C0OyIORcSYpM2SLqu7FoASmgS+R9IltvtsW9LVkkbrrgWghCbfgw9J\n2iRpm6QdE//NYOW9ABTQ6OfBI+JOSXdW3gVAYTwbBiRG4EBiBA4kRuBAYgQOJOaI8nen7H3Xijjn\njz9VfO7KNbuKz5Skp57sLz6zd/mx4jMl6eJle6rMXfuO7VXmXtl7qMrca/7qtuIzn7+sfAuSdN4t\nQ8VnDsV3dTQOu9NxnMGBxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgcSIzAgcQIHEiMwIHECBxIjMCB\nxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgcSIzAgcSq3FXV9iFJ/9Pg0HdK+mnxBeqZTfvOpl2l2bXv\nTNj13RHxK50OqhJ4U7aHI2KgtQWmaDbtO5t2lWbXvrNpVx6iA4kROJBY24EPtvz/n6rZtO9s2lWa\nXfvOml1b/R4cQF1tn8EBVNRa4Lavsf2M7Z22b29rj05sr7D9Pdujtkdsb2h7pyZsd9nebvtbbe8y\nGduLbG+y/aOJr/Glbe80Gdu3Tfw5eMr2122f1fZOk2klcNtdkr4g6VpJqySts72qjV0aOCnp0xHx\nXkmXSPrEDN71dBskjba9RAP3SHo4In5D0vs1g3e2vUzSrZIGImK1pC5JN7S71eTaOoNfLGlnROyK\niBOS7pe0tqVdJhURByNi28Svj2n8D+CydreanO3lkq6TdF/bu0zG9tmSrpD0JUmKiBMR8bN2t+qo\nW1Kv7W5JfZIOtLzPpNoKfJmkvad9vE8zPBpJst0v6UJJ5d/wuazPSfqMpFNtL9LBOZIOSfrKxLcT\n99me3/ZSZxIR+yXdJWmPpIOSXoqIb7e71eTaCvyN3rh8Rj+db/ttkr4h6ZMRcbTtfc7E9hpJL0TE\n1rZ3aaBb0kWS7o2ICyW9ImkmPx/zdo0/0lwpaamk+bZvbHerybUV+D5JK077eLlm8EMd23M1HvfG\niNjc9j4dXC7petvPavxbn6tsf63dlc5on6R9EfF/j4g2aTz4meoDknZHxKGIGJO0WdJlLe80qbYC\n3yLpPNsrbfdo/ImKb7a0y6RsW+PfI45GxN1t79NJRHw2IpZHRL/Gv66PRcSMPMtExHOS9to+f+JT\nV0t6usWVOtkj6RLbfRN/Lq7WDH5SUBp/iPSWi4iTtm+W9IjGn4n8ckSMtLFLA5dL+pikHbafnPjc\nHRHxUIs7ZXKLpI0Tf9HvknRTy/ucUUQM2d4kaZvGX13Zrhl+VRtXsgGJcSUbkBiBA4kROJAYgQOJ\nETiQGIEDiRE4kBiBA4n9L/yuy9QpJgQiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x8d4a6d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In python the input data starts at A55A980\n",
      "In C++ the input vector starts at address 000000000A87D5A0\n",
      "\n",
      "In python the reshaped data starts at A6B2040\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD8CAYAAABaQGkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAACnZJREFUeJzt3duvnXMex/HPZ3ar1RpxGDfaZpBg\nphGG7DgmLlSCIdzMBQnJuOnNOEYizI1/QIQLkTQON4SLciEiauJwMTcdWwlqIw2GOkSHjAoZdfjM\nxd6TlNG9nnY/P89+vt6vRNK9PZZPlvX2rLX26lMnEYCafjX0AADtEDhQGIEDhRE4UBiBA4UROFAY\ngQOFEThQGIEDhS1rcaMHeUVWanXvt3vCyV/1fputvPXKqqEnlPZLfyz8R19qT772pOPc4qOqh/qI\nnOENvd/ulg9f7v02W7ng6D8MPaG0X/pjYWue0e58NjFwnqIDhRE4UBiBA4UROFAYgQOFEThQWKfA\nbV9o+03bO2zf0noUgH5MDNz2lKS7JV0kab2kK2yvbz0MwOJ1OYOfLmlHkreT7JH0iKTL2s4C0Icu\nga+R9P5eX++c/94P2N5oe8b2zDf6uq99ABahS+A/9XG4//t8a5JNSaaTTC/XisUvA7BoXQLfKWnd\nXl+vlfRhmzkA+tQl8BckHW/7WNsHSbpc0uNtZwHow8TfLprkW9vXSNoiaUrS/Um2N18GYNE6/X7w\nJE9KerLxFgA945NsQGEEDhRG4EBhBA4URuBAYU2uqjo2Y7pAYquLDY7pPpDGt3conMGBwggcKIzA\ngcIIHCiMwIHCCBwojMCBwggcKIzAgcIIHCiMwIHCCBwojMCBwggcKIzAgcIIHCiMwIHCCBwojMCB\nwggcKIzAgcJGdVXVVlfSbHWl0hbGdh+M6eqnY3ocnH7BV52O4wwOFEbgQGEEDhRG4EBhBA4URuBA\nYRMDt73O9nO2Z21vt339zzEMwOJ1+Tn4t5JuSrLN9q8lvWj7b0leb7wNwCJNPIMn+SjJtvlffyFp\nVtKa1sMALN5+vQa3fYykUyVtbTEGQL86f1TV9iGSHpV0Q5LdP/H3N0raKEkrtaq3gQAOXKczuO3l\nmov7oSSP/dQxSTYlmU4yvVwr+twI4AB1eRfdku6TNJvkjvaTAPSlyxn8HElXSTrP9svzf/2x8S4A\nPZj4GjzJ3yX5Z9gCoGd8kg0ojMCBwggcKIzAgcIIHChsVBddHNNF8VrhPpgzpvuhxYUn38qnnY7j\nDA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEE\nDhRG4EBhBA4URuBAYQQOFOYkvd/o9Ckr848t63q/3VZaXPVybMZ0lVKJ/2Zb84x257OJf2YgZ3Cg\nMAIHCiNwoDACBwojcKAwAgcKI3CgsM6B256y/ZLtJ1oOAtCf/TmDXy9pttUQAP3rFLjttZIulnRv\n2zkA+tT1DH6npJslfb+vA2xvtD1je2bXp9/1Mg7A4kwM3PYlkj5J8uJCxyXZlGQ6yfRRR071NhDA\ngetyBj9H0qW235X0iKTzbD/YdBWAXkwMPMmtSdYmOUbS5ZKeTXJl82UAFo2fgwOFLdufg5M8L+n5\nJksA9I4zOFAYgQOFEThQGIEDhRE4UNh+vYuO7rhKaVst7t+x3QddcAYHCiNwoDACBwojcKAwAgcK\nI3CgMAIHCiNwoDACBwojcKAwAgcKI3CgMAIHCiNwoDACBwojcKAwAgcKI3CgMAIHCiNwoDACBwpr\nclXVt15ZNaorVHKFznZXgR3T/TCmK+GefsFXnY7jDA4URuBAYQQOFEbgQGEEDhRG4EBhnQK3fZjt\nzbbfsD1r+6zWwwAsXtefg98l6akkf7J9kKRVDTcB6MnEwG0fKulcSX+WpCR7JO1pOwtAH7o8RT9O\n0i5JD9h+yfa9tlc33gWgB10CXybpNEn3JDlV0peSbvnxQbY32p6xPfONvu55JoAD0SXwnZJ2Jtk6\n//VmzQX/A0k2JZlOMr1cK/rcCOAATQw8yceS3rd94vy3Nkh6vekqAL3o+i76tZIemn8H/W1JV7eb\nBKAvnQJP8rKk6cZbAPSMT7IBhRE4UBiBA4UROFAYgQOFEThQWJOrqp5w8lfasmU8Vyptcbtju0rp\nmK5+iu44gwOFEThQGIEDhRE4UBiBA4UROFAYgQOFEThQGIEDhRE4UBiBA4UROFAYgQOFEThQGIED\nhRE4UBiBA4UROFAYgQOFEThQWJOLLrbS6kKGLbS6iOGY7gOJizlKbe6Dt/Jpp+M4gwOFEThQGIED\nhRE4UBiBA4UROFAYgQOFdQrc9o22t9t+zfbDtle2HgZg8SYGbnuNpOskTSc5SdKUpMtbDwOweF2f\noi+TdLDtZZJWSfqw3SQAfZkYeJIPJN0u6T1JH0n6PMnTPz7O9kbbM7Zndn36Xf9LAey3Lk/RD5d0\nmaRjJR0tabXtK398XJJNSaaTTB915FT/SwHsty5P0c+X9E6SXUm+kfSYpLPbzgLQhy6BvyfpTNur\nbFvSBkmzbWcB6EOX1+BbJW2WtE3Sq/P/zKbGuwD0oNPvB09ym6TbGm8B0DM+yQYURuBAYQQOFEbg\nQGEEDhTmJL3f6KE+Imd4Q++3O7YriqKdFlcqbfX4arF1a57R7nzmScdxBgcKI3CgMAIHCiNwoDAC\nBwojcKAwAgcKI3CgMAIHCiNwoDACBwojcKAwAgcKI3CgMAIHCiNwoDACBwojcKAwAgcKI3CgMAIH\nCmtyVVXbuyT9s8Ohv5H0r94HtDOmvWPaKo1r71LY+tskR006qEngXdmeSTI92ID9NKa9Y9oqjWvv\nmLbyFB0ojMCBwoYOfNPA//79Naa9Y9oqjWvvaLYO+hocQFtDn8EBNDRY4LYvtP2m7R22bxlqxyS2\n19l+zvas7e22rx96Uxe2p2y/ZPuJobcsxPZhtjfbfmP+Pj5r6E0LsX3j/OPgNdsP21459KaFDBK4\n7SlJd0u6SNJ6SVfYXj/Elg6+lXRTkt9LOlPSX5bw1r1dL2l26BEd3CXpqSS/k3SKlvBm22skXSdp\nOslJkqYkXT7sqoUNdQY/XdKOJG8n2SPpEUmXDbRlQUk+SrJt/tdfaO4BuGbYVQuzvVbSxZLuHXrL\nQmwfKulcSfdJUpI9Sf497KqJlkk62PYySaskfTjwngUNFfgaSe/v9fVOLfFoJMn2MZJOlbR12CUT\n3SnpZknfDz1kguMk7ZL0wPzLiXttrx561L4k+UDS7ZLek/SRpM+TPD3sqoUNFfhP/cHlS/rtfNuH\nSHpU0g1Jdg+9Z19sXyLpkyQvDr2lg2WSTpN0T5JTJX0paSm/H3O45p5pHivpaEmrbV857KqFDRX4\nTknr9vp6rZbwUx3byzUX90NJHht6zwTnSLrU9ruae+lznu0Hh520Tzsl7Uzyv2dEmzUX/FJ1vqR3\nkuxK8o2kxySdPfCmBQ0V+AuSjrd9rO2DNPdGxeMDbVmQbWvuNeJskjuG3jNJkluTrE1yjObu12eT\nLMmzTJKPJb1v+8T5b22Q9PqAkyZ5T9KZtlfNPy42aAm/KSjNPUX62SX51vY1krZo7p3I+5NsH2JL\nB+dIukrSq7Zfnv/eX5M8OeCmSq6V9ND8/+jflnT1wHv2KclW25slbdPcT1de0hL/VBufZAMK45Ns\nQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhT2X771bWTJmQ2jAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x9ab35f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy\n",
    "import my_analysis_wrapper\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data = numpy.random.rand(10, 10, 10, 10)\n",
    "\n",
    "plt.imshow(data[:,:,5,5])\n",
    "plt.show()\n",
    "\n",
    "output = my_analysis_wrapper.process_with_vectors(data)\n",
    "plt.imshow(output[:,:,5,5])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Great, it worked\n",
    "\n",
    "## But we did copy our data  - several times\n",
    "\n",
    " - We copied the input once in Python to flatten it in the right order\n",
    " - Cython copied it again, because the pointer to the memory in the `vector` is different from the pointer to the data in the Numpy array.\n",
    " - Similarly, converting the output `vector` back to a Numpy array would involve a copy\n",
    " \n",
    "\n",
    "## Do we care?\n",
    "\n",
    "It depends on:\n",
    "\n",
    " - Is the processing time per-voxel much greater than the data copying time?\n",
    "   - If so, copying will not add significant overhead\n",
    " - Might the data be comparable in size to system memory?\n",
    "   - If so, copying may result in swapping and significant slowness\n",
    "   \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solution with less copying\n",
    "\n",
    "We can't use a `vector`, it needs to be free to manage its own memory, not use an existing fixed buffer.\n",
    "\n",
    "Instead pass a pure C array:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "    void process_array(float *input, int nx, int ny, int nz, int nt)\n",
    "    {\n",
    "        cerr << \"In C++ the input array starts at address \" << input << std::endl;\n",
    "\n",
    "        NEWIMAGE::volume4D<float> invol(nx, ny, nz, nt, input);\n",
    "\n",
    "        process_volume(invol);\n",
    "        \n",
    "        // Volume data buffer is modified directly, so provided it was not copied\n",
    "        // we should be able to see the output directly in Python   \n",
    "    }\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " - Note that we cannot check the size of the `input` buffer! It had better be correct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# my_analysis_wrapper.pyx\n",
    "\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "\n",
    "from libcpp.vector cimport vector\n",
    "\n",
    "cdef extern from \"my_analysis.h\":\n",
    "    void process_array(float *, int, int, int, int)\n",
    "\n",
    "def process_c(np.ndarray[np.float32_t, ndim=1] input,\n",
    "              nx, ny, nz, nt):\n",
    "    process_array(&input[0], nx, ny, nz, nt)\n",
    "\n",
    "def process_with_arrays(data):\n",
    "    # Save the dimension of the data because we're going to flatten it to 1D array\n",
    "    nx, ny, nz, nt = data.shape\n",
    "\n",
    "    # Convert data to 1D in Column-major (Fortran) order\n",
    "    data = data.flatten(order='F').astype(np.float32)\n",
    "\n",
    "    print(\"In python the data starts at %X\" % data.__array_interface__['data'][0])\n",
    "\n",
    "    process_c(data, nx, ny, nz, nt)\n",
    "\n",
    "    data = np.reshape(data, [nx, ny, nz, nt], order='F')\n",
    "    print(\"In python the reshaped data starts at %X\" % data.__array_interface__['data'][0])\n",
    "    \n",
    "    return data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD8CAYAAABaQGkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAADKRJREFUeJzt3W1sXvV5x/HfL3Yc7DQkTbulysPq\nwChrSFeBPMaDxgR0E5SMrGonkYpqY9WiTQVS2q2ibBKa9mJvECqqKjqPPkhrVMbSjFUdg9LSFqFp\nVpyEKhiXNktYHoGUpCTQ0Djk2gt7UsqI72P8/3PsS9+PhBSbw8UlJ9+c2/d9fG5HhADkNKftBQDU\nQ+BAYgQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJNZdY2jPnN7o7VpQfO6C839RfKYkdelU8Zkv7l9U\nfKYkvTavylh1v1rnisbXelxl7pIlR4rPPPDKwuIzJanr5fLn0RPHDuvk8Vc6fnGrBN7btUCXLv5I\n8blXPrCr+ExJWjDn1eIzN/71muIzJemllV1V5i5+ZqzK3JfePbfK3L/c8M/FZ9659Q+Kz5SkBU/0\nFZ/5kwfubnQcD9GBxAgcSIzAgcQIHEiMwIHECBxIrFHgtq+x/YztnbZvr70UgDI6Bm67S9IXJF0r\naZWkdbZX1V4MwPQ1OYNfLGlnROyKiBOS7pe0tu5aAEpoEvgySXtP+3jfxOd+ie31todtD584dbzU\nfgCmoUngb3S96/+7cDkiBiNiICIGeub0Tn8zANPWJPB9klac9vFySQfqrAOgpCaBb5F0nu2Vtnsk\n3SDpm3XXAlBCx58mi4iTtm+W9IikLklfjoiR6psBmLZGPy4aEQ9JeqjyLgAK40o2IDECBxIjcCAx\nAgcSI3AgsSo3XdTcbsWSdxQf+/TL5W+OKEkjX1xdfObLq+vcTfS3r9tRZe6WB99XZe7Y/Dp3a733\njvI39Tx/68HiMyVpbGn5r8Hunze7EzBncCAxAgcSI3AgMQIHEiNwIDECBxIjcCAxAgcSI3AgMQIH\nEiNwIDECBxIjcCAxAgcSI3AgMQIHEiNwIDECBxIjcCAxAgcSI3AgsTp3VZWk7vJ/d7zwh33FZ0rS\ndY/8oPjMB5/9zeIzJenx/7qgytwf3nx3lbm/9Y+fqjJ3/v7jxWf++3/WedPccx+7qfjMX/xNs+M4\ngwOJETiQGIEDiRE4kBiBA4kROJBYx8Btr7D9Pdujtkdsb3grFgMwfU1eBz8p6dMRsc32AklbbT8a\nEU9X3g3ANHU8g0fEwYjYNvHrY5JGJS2rvRiA6ZvS9+C2+yVdKGmoxjIAymocuO23SfqGpE9GxNE3\n+PfrbQ/bHj5x8ucldwTwJjUK3PZcjce9MSI2v9ExETEYEQMRMdDTXeeacQBT0+RZdEv6kqTRiKjz\nEwkAqmhyBr9c0sckXWX7yYl/Plh5LwAFdHyZLCKekOS3YBcAhXElG5AYgQOJETiQGIEDiRE4kFiV\nmy6+urhLO9ctLD73u+sGi8+UpFt3f7j4zL94z+PFZ0rS379Y5xXK9z36iSpz+4fGqsz98Z/OKz7z\nypG1xWdK0rJ/mVt85qEjzV7Y4gwOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGIEDiRG4EBiBA4k\nRuBAYgQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGIEDiRW5a6qPfPHtPyiA8Xn3nPoiuIzJWnX\nv51bfOaDG6P4TEla+rtVfst01vryv1+StHvdO6vMXflP5d8u76wf1nlf+1NHthWfOedks105gwOJ\nETiQGIEDiRE4kBiBA4kROJAYgQOJNQ7cdpft7ba/VXMhAOVM5Qy+QdJorUUAlNcocNvLJV0n6b66\n6wAoqekZ/HOSPiPp1JkOsL3e9rDt4bGf1bnkD8DUdAzc9hpJL0TE1smOi4jBiBiIiIG5i/qKLQjg\nzWtyBr9c0vW2n5V0v6SrbH+t6lYAiugYeER8NiKWR0S/pBskPRYRN1bfDMC08To4kNiUfrg4Ir4v\n6ftVNgFQHGdwIDECBxIjcCAxAgcSI3AgsSq36Bw72aUDhxcWn3vBoueKz5SkFf+6v/jMf9iyufhM\nSfrgtj+rMvfs639aZe7cO5ZVmfvcnx8rPrP30V8vPlOSeo6Vv8Pua//xRKPjOIMDiRE4kBiBA4kR\nOJAYgQOJETiQGIEDiRE4kBiBA4kROJAYgQOJETiQGIEDiRE4kBiBA4kROJAYgQOJETiQGIEDiRE4\nkBiBA4k5ovwdHxfOWxKXveujxec+//sris+UpMO/c6L4zPO+OFZ8piS9uHp+lbmrPj5SZe7zlx6t\nMvfjP95dfObQsXOLz5Sk31tY/mu7Ye1/6yc7jrvTcZzBgcQIHEiMwIHECBxIjMCBxAgcSKxR4LYX\n2d5k+0e2R21fWnsxANPX9N1F75H0cER8xHaPpL6KOwEopGPgts+WdIWkP5GkiDghqfyVIQCKa/IQ\n/RxJhyR9xfZ22/fZrnM5FYCimgTeLekiSfdGxIWSXpF0++sPsr3e9rDt4ROvHS+8JoA3o0ng+yTt\ni4ihiY83aTz4XxIRgxExEBEDPV29JXcE8CZ1DDwinpO01/b5E5+6WtLTVbcCUETTZ9FvkbRx4hn0\nXZJuqrcSgFIaBR4RT0oaqLwLgMK4kg1IjMCBxAgcSIzAgcQIHEiMwIHEmr4OPjX9UgyeKj72Vz9U\n5/qaJd85u/jMXXctKj5Tkvr/7kiVuUP9F1SZu7J7S5W5f7tjTfGZffPq3An34WffW3zmvuODjY7j\nDA4kRuBAYgQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGIE\nDiRG4EBiBA4kRuBAYlVuuvhr8w7r8+c8UHxu/0hf8ZmS9NWjS4vPvH/9tcVnStLO27uqzP3oBY9X\nmfvwM1dUmfvq3ig+c9sffb74TEn60PuvKT7zwJFmN4jkDA4kRuBAYgQOJEbgQGIEDiRG4EBiBA4k\n1ihw27fZHrH9lO2v2z6r9mIApq9j4LaXSbpV0kBErJbUJemG2osBmL6mD9G7JfXa7pbUJ+lAvZUA\nlNIx8IjYL+kuSXskHZT0UkR8+/XH2V5ve9j28OHD5d8bHMDUNXmI/nZJayWtlLRU0nzbN77+uIgY\njIiBiBhYvJjn7oCZoEmJH5C0OyIORcSYpM2SLqu7FoASmgS+R9IltvtsW9LVkkbrrgWghCbfgw9J\n2iRpm6QdE//NYOW9ABTQ6OfBI+JOSXdW3gVAYTwbBiRG4EBiBA4kRuBAYgQOJOaI8nen7H3Xijjn\njz9VfO7KNbuKz5Skp57sLz6zd/mx4jMl6eJle6rMXfuO7VXmXtl7qMrca/7qtuIzn7+sfAuSdN4t\nQ8VnDsV3dTQOu9NxnMGBxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgcSIzAgcQIHEiMwIHECBxIjMCB\nxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgcSIzAgcSq3FXV9iFJ/9Pg0HdK+mnxBeqZTfvOpl2l2bXv\nTNj13RHxK50OqhJ4U7aHI2KgtQWmaDbtO5t2lWbXvrNpVx6iA4kROJBY24EPtvz/n6rZtO9s2lWa\nXfvOml1b/R4cQF1tn8EBVNRa4Lavsf2M7Z22b29rj05sr7D9Pdujtkdsb2h7pyZsd9nebvtbbe8y\nGduLbG+y/aOJr/Glbe80Gdu3Tfw5eMr2122f1fZOk2klcNtdkr4g6VpJqySts72qjV0aOCnp0xHx\nXkmXSPrEDN71dBskjba9RAP3SHo4In5D0vs1g3e2vUzSrZIGImK1pC5JN7S71eTaOoNfLGlnROyK\niBOS7pe0tqVdJhURByNi28Svj2n8D+CydreanO3lkq6TdF/bu0zG9tmSrpD0JUmKiBMR8bN2t+qo\nW1Kv7W5JfZIOtLzPpNoKfJmkvad9vE8zPBpJst0v6UJJ5d/wuazPSfqMpFNtL9LBOZIOSfrKxLcT\n99me3/ZSZxIR+yXdJWmPpIOSXoqIb7e71eTaCvyN3rh8Rj+db/ttkr4h6ZMRcbTtfc7E9hpJL0TE\n1rZ3aaBb0kWS7o2ICyW9ImkmPx/zdo0/0lwpaamk+bZvbHerybUV+D5JK077eLlm8EMd23M1HvfG\niNjc9j4dXC7petvPavxbn6tsf63dlc5on6R9EfF/j4g2aTz4meoDknZHxKGIGJO0WdJlLe80qbYC\n3yLpPNsrbfdo/ImKb7a0y6RsW+PfI45GxN1t79NJRHw2IpZHRL/Gv66PRcSMPMtExHOS9to+f+JT\nV0t6usWVOtkj6RLbfRN/Lq7WDH5SUBp/iPSWi4iTtm+W9IjGn4n8ckSMtLFLA5dL+pikHbafnPjc\nHRHxUIs7ZXKLpI0Tf9HvknRTy/ucUUQM2d4kaZvGX13Zrhl+VRtXsgGJcSUbkBiBA4kROJAYgQOJ\nETiQGIEDiRE4kBiBA4n9L/yuy9QpJgQiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa4908d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In python the data starts at A870040\n",
      "In C++ the input array starts at address 000000000A870040\n",
      "\n",
      "In python the reshaped data starts at A870040\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD8CAYAAABaQGkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAACnZJREFUeJzt3duvnXMex/HPZ3ar1RpxGDfaZpBg\nphGG7DgmLlSCIdzMBQnJuOnNOEYizI1/QIQLkTQON4SLciEiauJwMTcdWwlqIw2GOkSHjAoZdfjM\nxd6TlNG9nnY/P89+vt6vRNK9PZZPlvX2rLX26lMnEYCafjX0AADtEDhQGIEDhRE4UBiBA4UROFAY\ngQOFEThQGIEDhS1rcaMHeUVWanXvt3vCyV/1fputvPXKqqEnlPZLfyz8R19qT772pOPc4qOqh/qI\nnOENvd/ulg9f7v02W7ng6D8MPaG0X/pjYWue0e58NjFwnqIDhRE4UBiBA4UROFAYgQOFEThQWKfA\nbV9o+03bO2zf0noUgH5MDNz2lKS7JV0kab2kK2yvbz0MwOJ1OYOfLmlHkreT7JH0iKTL2s4C0Icu\nga+R9P5eX++c/94P2N5oe8b2zDf6uq99ABahS+A/9XG4//t8a5JNSaaTTC/XisUvA7BoXQLfKWnd\nXl+vlfRhmzkA+tQl8BckHW/7WNsHSbpc0uNtZwHow8TfLprkW9vXSNoiaUrS/Um2N18GYNE6/X7w\nJE9KerLxFgA945NsQGEEDhRG4EBhBA4URuBAYU2uqjo2Y7pAYquLDY7pPpDGt3conMGBwggcKIzA\ngcIIHCiMwIHCCBwojMCBwggcKIzAgcIIHCiMwIHCCBwojMCBwggcKIzAgcIIHCiMwIHCCBwojMCB\nwggcKIzAgcJGdVXVVlfSbHWl0hbGdh+M6eqnY3ocnH7BV52O4wwOFEbgQGEEDhRG4EBhBA4URuBA\nYRMDt73O9nO2Z21vt339zzEMwOJ1+Tn4t5JuSrLN9q8lvWj7b0leb7wNwCJNPIMn+SjJtvlffyFp\nVtKa1sMALN5+vQa3fYykUyVtbTEGQL86f1TV9iGSHpV0Q5LdP/H3N0raKEkrtaq3gQAOXKczuO3l\nmov7oSSP/dQxSTYlmU4yvVwr+twI4AB1eRfdku6TNJvkjvaTAPSlyxn8HElXSTrP9svzf/2x8S4A\nPZj4GjzJ3yX5Z9gCoGd8kg0ojMCBwggcKIzAgcIIHChsVBddHNNF8VrhPpgzpvuhxYUn38qnnY7j\nDA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEE\nDhRG4EBhBA4URuBAYQQOFOYkvd/o9Ckr848t63q/3VZaXPVybMZ0lVKJ/2Zb84x257OJf2YgZ3Cg\nMAIHCiNwoDACBwojcKAwAgcKI3CgsM6B256y/ZLtJ1oOAtCf/TmDXy9pttUQAP3rFLjttZIulnRv\n2zkA+tT1DH6npJslfb+vA2xvtD1je2bXp9/1Mg7A4kwM3PYlkj5J8uJCxyXZlGQ6yfRRR071NhDA\ngetyBj9H0qW235X0iKTzbD/YdBWAXkwMPMmtSdYmOUbS5ZKeTXJl82UAFo2fgwOFLdufg5M8L+n5\nJksA9I4zOFAYgQOFEThQGIEDhRE4UNh+vYuO7rhKaVst7t+x3QddcAYHCiNwoDACBwojcKAwAgcK\nI3CgMAIHCiNwoDACBwojcKAwAgcKI3CgMAIHCiNwoDACBwojcKAwAgcKI3CgMAIHCiNwoDACBwpr\nclXVt15ZNaorVHKFznZXgR3T/TCmK+GefsFXnY7jDA4URuBAYQQOFEbgQGEEDhRG4EBhnQK3fZjt\nzbbfsD1r+6zWwwAsXtefg98l6akkf7J9kKRVDTcB6MnEwG0fKulcSX+WpCR7JO1pOwtAH7o8RT9O\n0i5JD9h+yfa9tlc33gWgB10CXybpNEn3JDlV0peSbvnxQbY32p6xPfONvu55JoAD0SXwnZJ2Jtk6\n//VmzQX/A0k2JZlOMr1cK/rcCOAATQw8yceS3rd94vy3Nkh6vekqAL3o+i76tZIemn8H/W1JV7eb\nBKAvnQJP8rKk6cZbAPSMT7IBhRE4UBiBA4UROFAYgQOFEThQWJOrqp5w8lfasmU8Vyptcbtju0rp\nmK5+iu44gwOFEThQGIEDhRE4UBiBA4UROFAYgQOFEThQGIEDhRE4UBiBA4UROFAYgQOFEThQGIED\nhRE4UBiBA4UROFAYgQOFEThQWJOLLrbS6kKGLbS6iOGY7gOJizlKbe6Dt/Jpp+M4gwOFEThQGIED\nhRE4UBiBA4UROFAYgQOFdQrc9o22t9t+zfbDtle2HgZg8SYGbnuNpOskTSc5SdKUpMtbDwOweF2f\noi+TdLDtZZJWSfqw3SQAfZkYeJIPJN0u6T1JH0n6PMnTPz7O9kbbM7Zndn36Xf9LAey3Lk/RD5d0\nmaRjJR0tabXtK398XJJNSaaTTB915FT/SwHsty5P0c+X9E6SXUm+kfSYpLPbzgLQhy6BvyfpTNur\nbFvSBkmzbWcB6EOX1+BbJW2WtE3Sq/P/zKbGuwD0oNPvB09ym6TbGm8B0DM+yQYURuBAYQQOFEbg\nQGEEDhTmJL3f6KE+Imd4Q++3O7YriqKdFlcqbfX4arF1a57R7nzmScdxBgcKI3CgMAIHCiNwoDAC\nBwojcKAwAgcKI3CgMAIHCiNwoDACBwojcKAwAgcKI3CgMAIHCiNwoDACBwojcKAwAgcKI3CgMAIH\nCmtyVVXbuyT9s8Ohv5H0r94HtDOmvWPaKo1r71LY+tskR006qEngXdmeSTI92ID9NKa9Y9oqjWvv\nmLbyFB0ojMCBwoYOfNPA//79Naa9Y9oqjWvvaLYO+hocQFtDn8EBNDRY4LYvtP2m7R22bxlqxyS2\n19l+zvas7e22rx96Uxe2p2y/ZPuJobcsxPZhtjfbfmP+Pj5r6E0LsX3j/OPgNdsP21459KaFDBK4\n7SlJd0u6SNJ6SVfYXj/Elg6+lXRTkt9LOlPSX5bw1r1dL2l26BEd3CXpqSS/k3SKlvBm22skXSdp\nOslJkqYkXT7sqoUNdQY/XdKOJG8n2SPpEUmXDbRlQUk+SrJt/tdfaO4BuGbYVQuzvVbSxZLuHXrL\nQmwfKulcSfdJUpI9Sf497KqJlkk62PYySaskfTjwngUNFfgaSe/v9fVOLfFoJMn2MZJOlbR12CUT\n3SnpZknfDz1kguMk7ZL0wPzLiXttrx561L4k+UDS7ZLek/SRpM+TPD3sqoUNFfhP/cHlS/rtfNuH\nSHpU0g1Jdg+9Z19sXyLpkyQvDr2lg2WSTpN0T5JTJX0paSm/H3O45p5pHivpaEmrbV857KqFDRX4\nTknr9vp6rZbwUx3byzUX90NJHht6zwTnSLrU9ruae+lznu0Hh520Tzsl7Uzyv2dEmzUX/FJ1vqR3\nkuxK8o2kxySdPfCmBQ0V+AuSjrd9rO2DNPdGxeMDbVmQbWvuNeJskjuG3jNJkluTrE1yjObu12eT\nLMmzTJKPJb1v+8T5b22Q9PqAkyZ5T9KZtlfNPy42aAm/KSjNPUX62SX51vY1krZo7p3I+5NsH2JL\nB+dIukrSq7Zfnv/eX5M8OeCmSq6V9ND8/+jflnT1wHv2KclW25slbdPcT1de0hL/VBufZAMK45Ns\nQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEEDhT2X771bWTJmQ2jAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xaa4b128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy\n",
    "import my_analysis_wrapper\n",
    "\n",
    "plt.imshow(data[:,:,5,5])\n",
    "plt.show()\n",
    "\n",
    "output = my_analysis_wrapper.process_with_arrays(data)\n",
    "plt.imshow(output[:,:,5,5])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " - We copied our input data once when we flattened it into Fortran order\n",
    " - C++ code operated directly on that buffer\n",
    " - Output data was not copied when reshaped "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary\n",
    "\n",
    " - Easy-ish recipe for passing Numpy arrays to C++ either as a `std::vector` or as a `float *` array.\n",
    " - Can construct `NEWIMAGE::volume<float>` or other complex containers from within C++\n",
    " - Easy modification to instead use `double` array\n",
    " - Can pass Python strings to `C++ std::string` and other C++ containers in a similar way\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Alternatives (Briefly!)\n",
    "\n",
    "## Why?\n",
    "\n",
    "![Compatibility](cython_compat1.png)\n",
    "\n",
    "## Can we assume that our newly compiled Cython/C++ code will link correctly with `libnewimage.a`?\n",
    "\n",
    " - Often, yes, but in general, no\n",
    " - It depends on the compiler used for each - ideally they need to match\n",
    " - The compiler of your Cython extension is **fixed** by the version of Python you are using\n",
    " - Might need to recompile your dependency libraries with this compiler\n",
    " - If you can't do this (e.g. commercial binary) you may be **stuck**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Three common problem scenarios\n",
    "\n",
    " - On Mac, need to use the same C++ standard library (either `libc++` or `libstdc++`)\n",
    " - On Python 2, C++ compiler will be very old (may not support all of C++11)\n",
    " - On Windows, no two versions of VC++ are binary compatible\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Alternative approach where this is a problem\n",
    "\n",
    " - Make your code a shared library with a *Pure C* API\n",
    " - Use `ctypes`\n",
    " \n",
    "## `ctypes`\n",
    "\n",
    " - Part of Python standard library\n",
    " - Allows you to call library functions from 'C' shared library (not C++)\n",
    " - **Pure 'C' libraries are (generally) binary compatible on a given platform**\n",
    " - We have to load the library manually\n",
    " - We have to tell Python about the input and return types\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pure 'C' API for our processing function\n",
    "\n",
    "    // my_analysis_purec.h\n",
    "\n",
    "    #ifdef __cplusplus\n",
    "    extern \"C\" {\n",
    "    #endif\n",
    "    \n",
    "    void process_array(float *input, int nx, int ny, int nz, int nt);\n",
    "    \n",
    "    #ifdef __cplusplus\n",
    "    }\n",
    "    #endif\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " - Note need to use `extern \"C\" { }` if we may want to include this header from C++\n",
    " - **On Windows, additional code is required to make the shared library (DLL) link correctly!**\n",
    " - Note that the implementation *can* use C++"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import ctypes import CDLL, c_int, c_char_p\n",
    "import numpy as np\n",
    "import numpy.ctypeslib\n",
    "\n",
    "def process_ctypes(data):\n",
    "    \n",
    "    clib = ctypes.cdll.LoadLibrary(\"libmy_analysis.so\")\n",
    "\n",
    "    # This is the data type of a 1-D Numpy array\n",
    "    c_float_arr = numpy.ctypeslib.ndpointer(dtype=np.float32, ndim=1, flags='CONTIGUOUS')\n",
    "\n",
    "    # This specifies the argument types for the 'process_array' function\n",
    "    # This is not actually required but enables ctypes to do some error checking\n",
    "    clib.process_array.argtypes = [c_float_arr, c_int, c_int, c_int, c_int]\n",
    "\n",
    "    # Put the Numpy data into row-major order and make sure it is contiguous in memory\n",
    "    item = np.ascontiguousarray(item.flatten(order='F'), dtype=np.float32)\n",
    "    \n",
    "    clib.process_carray(data, shape[0], shape[1], shape[2], shape[3])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison with Cython\n",
    "\n",
    "## Cython advantages\n",
    "\n",
    " - Python wrapper is probably a little quicker and cleaner to write\n",
    " - Don't need to produce a new pure-C API provided we have an entry point using C++ types\n",
    " - Potential for better error-checking\n",
    " - Might integrate well if you are already using Cython\n",
    " - No need to build a shared library\n",
    " \n",
    "## `ctypes` advantages\n",
    "\n",
    " - Part of the Python standard library\n",
    " - No additional compile step in `setup.py`\n",
    " - Binary compatibility - no need to be tied to a single (perhaps old) C++ compiler\n",
    " \n",
    "## Conclusion?\n",
    "\n",
    " - Use Cython when you can, `ctypes` if you have to\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Other alternatives (briefly for completeness)\n",
    "\n",
    "## Wrapper Generators (SWIG, shiboken, others)\n",
    "\n",
    " - Run a preprocessor on your C++ code to generate an 'automatic' Python wrapper\n",
    " - Generally need to write an 'interface specifier' for each C++ header to describe how it interfaces to Python\n",
    " - Great when you have a large, complex C++ API which needs to be consistently exposed to Python (e.g. wx/wxpython, QT/PyQT)\n",
    " - SWIG can support other languages as well as Python\n",
    " - Probably more work than Cython/ctypes if you have a single simple API\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Just Wrap the Command Line\n",
    "\n",
    " - Quick and dirty\n",
    " - Copies all data to/from filesystem\n",
    " - Need to go via command line API, create temp directories, etc\n",
    " - Don't overlook as a way of getting started - can move to other solution later\n",
    " \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD8CAYAAABaQGkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAADJxJREFUeJzt3V+s1/V9x/HXi3OAI+eUwopzCKxC\ntVps12BPOiuJWcVk/WOlSXdhV01nlpLp6p/WrLVNGrdkF72wziY6EkI16XC6hnJhq1X7x27r6kiP\nwLR4SnqGFU9FRQGhWMAD716cs4Q6Ob/v4Xw+fs955/lISDiHr2/fIefJ93d+53c+xxEhADnNaHsB\nAPUQOJAYgQOJETiQGIEDiRE4kBiBA4kROJAYgQOJddcY2jOvJ/oW9hWfO7/7UPGZkjT84oLiM6Pv\nePGZkjRjRp253UNHqsyde/6xKnMP7phdfOay5fuLz5Sk/32ifAuHdUhH44g7XVcl8L6FfVr9zcuK\nz/3E2waKz5SkL97x18Vn/vZP6/xj1Nd7uMrc0z8+VGXupd+qE82jH1xWfOa3vved4jMl6ROLLyw+\nc3P8sNF1PEQHEiNwIDECBxIjcCAxAgcSI3AgsUaB2/6Q7R22h2zfXHspAGV0DNx2l6Q7JX1Y0nJJ\nn7S9vPZiACavyR38/ZKGImJnRByVdJ+k1XXXAlBCk8AXSXr2hLeHx973e2yvsT1ge+Dw/jqvtgIw\nMU0Cf6PXu/6/o1gjYl1E9EdEf8+8nslvBmDSmgQ+LGnJCW8vlvRcnXUAlNQk8J9JOsf2UtuzJF0h\n6f66awEooeN3k0XEiO3PSnpYUpekuyJie/XNAExao28XjYgHJT1YeRcAhfFKNiAxAgcSI3AgMQIH\nEiNwIDHX+Png576nJ9be//bic796Xn/xmZJ09k87Hk45YT984H3FZ0rS0bl1TlX9o8fq/Jz4tz7x\ncpW5g1+cV3zmoge6is+UpNtuvaP4zKs/tluDT3Q+VZU7OJAYgQOJETiQGIEDiRE4kBiBA4kROJAY\ngQOJETiQGIEDiRE4kBiBA4kROJAYgQOJETiQGIEDiRE4kBiBA4kROJAYgQOJETiQWKOfTTZR+471\natO+8qeKzji7/EmtkjR0zeziMz+4dkvxmZL0zOXlTxOVpKf+fknni07BW7ceqzJ36X3lZ8758q7y\nQyXtOfaW4jNH4oVG13EHBxIjcCAxAgcSI3AgMQIHEiNwILGOgdteYvtR24O2t9u+4c1YDMDkNfk6\n+IikmyJii+23SHrc9vcj4qnKuwGYpI538IjYHRFbxn5/UNKgpEW1FwMweRP6HNz2WZJWSNpcYxkA\nZTUO3HafpG9LujEiDrzBn6+xPWB74Lf7DpfcEcApahS47ZkajfueiNj0RtdExLqI6I+I/tPm95Tc\nEcApavIsuiV9Q9JgRNxWfyUApTS5g6+UdJWkS2xvG/v1kcp7ASig45fJIuInkvwm7AKgMF7JBiRG\n4EBiBA4kRuBAYgQOJFbl0MUjx7u18zcLyg8eqXOAXw3/+W8XVJk78+NRZ+7LVcbq6U8trDL3tbnH\ni8+8av4zxWdK0o7DZxafeTiaHRDJHRxIjMCBxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgcSIzAgcQI\nHEiMwIHECBxIjMCBxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgcSKzKqaqzZozoj3v3FZ/7yD+cUXym\nJC355sziM5/8/D8XnylJ/bdcU2Vuz0t1fvzcqwvrnAL7zrUvFJ+5cfefFZ8pSX2ryu+69+gTja7j\nDg4kRuBAYgQOJEbgQGIEDiRG4EBiBA4k1jhw2122t9r+bs2FAJQzkTv4DZIGay0CoLxGgdteLOmj\nktbXXQdASU3v4LdL+oKkk/7UddtrbA/YHji870iR5QBMTsfAbV8m6cWIeHy86yJiXUT0R0R/z/zZ\nxRYEcOqa3MFXSrrc9q8k3SfpEtsbqm4FoIiOgUfElyJicUScJekKST+KiCurbwZg0vg6OJDYhL4f\nPCJ+LOnHVTYBUBx3cCAxAgcSI3AgMQIHEiNwILEqp6ou7D6gm8/4QfG5axb8e/GZkvT9P1lefObS\nBz5TfKYkzf/Y3ipzjz32tipzZ++tc1rryB/OLT5z7qXPF58pSa/dW+E04L3N0uUODiRG4EBiBA4k\nRuBAYgQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGIEDiRG4EBiBA4kRuBAYgQOJEbgQGIEDiRG\n4EBiBA4kVuVU1ZeP9WrD/vcVn/tfq88rPlOS7v6Pfy0+87FzlxWfKUnbttWZe8Yzx6vM/ad/vLPK\n3JuGri0+84VKJ8B2v6P8zGOzm13HHRxIjMCBxAgcSIzAgcQIHEiMwIHEGgVue57tjbZ/YXvQ9gdq\nLwZg8pp+Hfzrkh6KiL+wPUvSnIo7ASikY+C250q6WNJfSVJEHJV0tO5aAEpo8hB9maQ9ku62vdX2\netu9lfcCUECTwLslXSBpbUSskHRI0s2vv8j2GtsDtgcO7eMGD0wFTQIfljQcEZvH3t6o0eB/T0Ss\ni4j+iOjvnT+r5I4ATlHHwCPieUnP2j537F2rJD1VdSsARTR9Fv06SfeMPYO+U9LV9VYCUEqjwCNi\nm6T+yrsAKIxXsgGJETiQGIEDiRE4kBiBA4kROJCYI6L40NPOWBJn/+Xni8898J46L4FduXyo+Myf\n/rLO6ad9/9NTZe7CnxysMveVd9b5toVDZ5a/N824aF/xmZLU3XWs+MwdN96lV3+5u+MxsNzBgcQI\nHEiMwIHECBxIjMCBxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgc\nSIzAgcQIHEis6Q8fnJAFpx/Qpz/zUPG5G+748+IzJenAOeUPMuzdXudwxIPn1zl48tCKKh8K+t7F\nX6sy99pPX1d85ux/qXPoYhw/XnzmzpeaHeTIHRxIjMCBxAgcSIzAgcQIHEiMwIHECBxIrFHgtj9n\ne7vtn9u+13adL/ICKKpj4LYXSbpeUn9EvFtSl6Qrai8GYPKaPkTvlnSa7W5JcyQ9V28lAKV0DDwi\nfi3pVkm7JO2W9EpEPPL662yvsT1ge+A3e+u8nBLAxDR5iD5f0mpJSyWdKanX9pWvvy4i1kVEf0T0\n9/3BrPKbApiwJg/RL5X0dETsiYjXJG2SdFHdtQCU0CTwXZIutD3HtiWtkjRYdy0AJTT5HHyzpI2S\ntkh6cuy/WVd5LwAFNPom4Ii4RdItlXcBUBivZAMSI3AgMQIHEiNwIDECBxKrcpTmDB1XX9fh4nP3\nX3Sk+ExJOvLwsuIz1/3NHcVnStJNX7m2ytyZr3ZVmXvN3ddXmft36zcUn3n7u95bfKYknf/fI8Vn\nbv1Us5ncwYHECBxIjMCBxAgcSIzAgcQIHEiMwIHECBxIjMCBxAgcSIzAgcQIHEiMwIHECBxIjMCB\nxAgcSIzAgcQIHEiMwIHECBxIjMCBxBwR5YfaeyQ90+DSBZJeKr5APdNp3+m0qzS99p0Ku749Ik7v\ndFGVwJuyPRAR/a0tMEHTad/ptKs0vfadTrvyEB1IjMCBxNoOfF3L//+Jmk77Tqddpem177TZtdXP\nwQHU1fYdHEBFrQVu+0O2d9gesn1zW3t0YnuJ7UdtD9rebvuGtndqwnaX7a22v9v2LuOxPc/2Rtu/\nGPs7/kDbO43H9ufGPg5+bvte2z1t7zSeVgK33SXpTkkflrRc0idtL29jlwZGJN0UEe+SdKGkv53C\nu57oBkmDbS/RwNclPRQR50l6r6bwzrYXSbpeUn9EvFtSl6Qr2t1qfG3dwd8vaSgidkbEUUn3SVrd\n0i7jiojdEbFl7PcHNfoBuKjdrcZne7Gkj0pa3/Yu47E9V9LFkr4hSRFxNCL2t7tVR92STrPdLWmO\npOda3mdcbQW+SNKzJ7w9rCkejSTZPkvSCkmb292ko9slfUHS8bYX6WCZpD2S7h77dGK97d62lzqZ\niPi1pFsl7ZK0W9IrEfFIu1uNr63A/Qbvm9JP59vuk/RtSTdGxIG29zkZ25dJejEiHm97lwa6JV0g\naW1ErJB0SNJUfj5mvkYfaS6VdKakXttXtrvV+NoKfFjSkhPeXqwp/FDH9kyNxn1PRGxqe58OVkq6\n3PavNPqpzyW2N7S70kkNSxqOiP97RLRRo8FPVZdKejoi9kTEa5I2Sbqo5Z3G1VbgP5N0ju2ltmdp\n9ImK+1vaZVy2rdHPEQcj4ra29+kkIr4UEYsj4iyN/r3+KCKm5F0mIp6X9Kztc8fetUrSUy2u1Mku\nSRfanjP2cbFKU/hJQWn0IdKbLiJGbH9W0sMafSbyrojY3sYuDayUdJWkJ21vG3vflyPiwRZ3yuQ6\nSfeM/UO/U9LVLe9zUhGx2fZGSVs0+tWVrZrir2rjlWxAYrySDUiMwIHECBxIjMCBxAgcSIzAgcQI\nHEiMwIHEfgdwPctS4OQ+QAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xb047588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD8CAYAAABaQGkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAACn1JREFUeJzt3c+v3XMex/HXa25JlRFUN9pmSmLM\niIwhJx0/EguVYAibWZCQjE03gxKJMBv/gAgLkTR+bAiLshARNfFjMZtylRnqMhGM1o+4UxkVYkq9\nZnHvJGXae77t/X5873l7PhJJ73Ucr5zep++5555+OIkA1PSzoQcAaIfAgcIIHCiMwIHCCBwojMCB\nwggcKIzAgcIIHChsWYs7PfGEqaxbe0SLu27iH39f0ft9/vI3X/V+ny21eAykdo/DJP2etdj6tb7U\n3vzH427nFm9VHZ25PC9tXdv7/bZy8Um/7f0+t370Wu/32VKLx0Bq9zhM0u9Zi63b8pz25LOxgfMU\nHSiMwIHCCBwojMCBwggcKIzAgcI6BW77Ettv237H9m2tRwHox9jAbU9JulfSpZJOl3S17dNbDwOw\neF2u4OslvZPk3SR7JT0m6cq2swD0oUvgqyXt3O/jXfOf+x7bG21P256e3b2vr30AFqFL4Ad6O9z/\nvb81yeYkoySjVSunFr8MwKJ1CXyXpP3fWL5G0kdt5gDoU5fAX5Z0qu2TbR8p6SpJT7adBaAPY/+4\naJJvbV8vaaukKUkPJtnRfBmARev058GTPC3p6cZbAPSMd7IBhRE4UBiBA4UROFAYgQOFNTlVtZVJ\nOhhwkra2NEmPwyRtXX9xtxNguYIDhRE4UBiBA4UROFAYgQOFEThQGIEDhRE4UBiBA4UROFAYgQOF\nEThQGIEDhRE4UBiBA4UROFAYgQOFEThQGIEDhRE4UBiBA4VN1KmqrbQ4TbPV6aeTdPKn1G7vJP2e\nDYkrOFAYgQOFEThQGIEDhRE4UBiBA4WNDdz2Wtsv2J6xvcP2ph9jGIDF6/Jz8G8l3ZJku+2fS3rF\n9l+SvNl4G4BFGnsFT/Jxku3zv/5C0oyk1a2HAVi8Q/oe3PY6SWdJ2tZiDIB+dQ7c9jGSHpd0U5I9\nB/j7G21P256e3b2vz40ADlOnwG0fobm4H0nyxIFuk2RzklGS0aqVU31uBHCYuryKbkkPSJpJclf7\nSQD60uUKfr6kayVdaPu1+b9+33gXgB6M/TFZkr9K8o+wBUDPeCcbUBiBA4UROFAYgQOFEThQGIcu\nNtLqsMFWJm1vxQMSW+AKDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4URuBAYQQOFEbgQGEE\nDhRG4EBhBA4URuBAYQQOFEbgQGEEDhRG4EBhBA4UNlGnqrY6SbPFiaKTtLUlHodhT4DlCg4URuBA\nYQQOFEbgQGEEDhRG4EBhBA4U1jlw21O2X7X9VMtBAPpzKFfwTZJmWg0B0L9OgdteI+kySfe3nQOg\nT12v4HdLulXSdwe7ge2NtqdtT8/u3tfLOACLMzZw25dL+jTJKwvdLsnmJKMko1Urp3obCODwdbmC\nny/pCtvvS3pM0oW2H266CkAvxgae5PYka5Ksk3SVpOeTXNN8GYBF4+fgQGGH9OfBk7wo6cUmSwD0\njis4UBiBA4UROFAYgQOFEThQmJP0fqejM5fnpa1re79ftDNJp5S2MkknwG7Lc9qTzzzudlzBgcII\nHCiMwIHCCBwojMCBwggcKIzAgcIIHCiMwIHCCBwojMCBwggcKIzAgcIIHCiMwIHCCBwojMCBwggc\nKIzAgcIIHCiMwIHCDun/TTa0Vid/tjpNE5N1UmlFXMGBwggcKIzAgcIIHCiMwIHCCBworFPgto+z\nvcX2W7ZnbJ/behiAxev6c/B7JD2T5A+2j5S0ouEmAD0ZG7jtYyVdIOmPkpRkr6S9bWcB6EOXp+in\nSJqV9JDtV23fb/voxrsA9KBL4MsknS3pviRnSfpS0m0/vJHtjbanbU/P7t7X80wAh6NL4Lsk7Uqy\nbf7jLZoL/nuSbE4ySjJatXKqz40ADtPYwJN8Immn7dPmP7VB0ptNVwHoRddX0W+Q9Mj8K+jvSrqu\n3SQAfekUeJLXJI0abwHQM97JBhRG4EBhBA4URuBAYQQOFEbgQGFO0vudHusT8jtv6P1+Of2U00Rb\nmqSvr/UX79T03772uNtxBQcKI3CgMAIHCiNwoDACBwojcKAwAgcKI3CgMAIHCiNwoDACBwojcKAw\nAgcKI3CgMAIHCiNwoDACBwojcKAwAgcKI3CgsCaHLo7OXJ6Xtq7t/X5bHTjY4rC9Sdo6iX7qh09u\ny3Pak884dBH4KSNwoDACBwojcKAwAgcKI3CgMAIHCusUuO2bbe+w/YbtR20vbz0MwOKNDdz2akk3\nSholOUPSlKSrWg8DsHhdn6Ivk3SU7WWSVkj6qN0kAH0ZG3iSDyXdKekDSR9L+jzJsz+8ne2Ntqdt\nT8/u3tf/UgCHrMtT9OMlXSnpZEknSTra9jU/vF2SzUlGSUarVk71vxTAIevyFP0iSe8lmU3yjaQn\nJJ3XdhaAPnQJ/ANJ59heYduSNkiaaTsLQB+6fA++TdIWSdslvT7/z2xuvAtAD5Z1uVGSOyTd0XgL\ngJ7xTjagMAIHCiNwoDACBwojcKCwTq+iLxWtThRtcULnJG2dRD/1k3DXX/xVp9txBQcKI3CgMAIH\nCiNwoDACBwojcKAwAgcKI3CgMAIHCiNwoDACBwojcKAwAgcKI3CgMAIHCiNwoDACBwojcKAwAgcK\nI3CgMAIHCnOS/u/UnpX0zw43PVHSv3of0M4k7Z2krdJk7V0KW3+RZNW4GzUJvCvb00lGgw04RJO0\nd5K2SpO1d5K28hQdKIzAgcKGDnzzwP/+QzVJeydpqzRZeydm66DfgwNoa+grOICGBgvc9iW237b9\nju3bhtoxju21tl+wPWN7h+1NQ2/qwvaU7VdtPzX0loXYPs72FttvzT/G5w69aSG2b57/OnjD9qO2\nlw+9aSGDBG57StK9ki6VdLqkq22fPsSWDr6VdEuSX0s6R9KflvDW/W2SNDP0iA7ukfRMkl9JOlNL\neLPt1ZJulDRKcoakKUlXDbtqYUNdwddLeifJu0n2SnpM0pUDbVlQko+TbJ//9Rea+wJcPeyqhdle\nI+kySfcPvWUhto+VdIGkByQpyd4k/x521VjLJB1le5mkFZI+GnjPgoYKfLWknft9vEtLPBpJsr1O\n0lmStg27ZKy7Jd0q6buhh4xxiqRZSQ/Nfztxv+2jhx51MEk+lHSnpA8kfSzp8yTPDrtqYUMF7gN8\nbkm/nG/7GEmPS7opyZ6h9xyM7cslfZrklaG3dLBM0tmS7ktylqQvJS3l12OO19wzzZMlnSTpaNvX\nDLtqYUMFvkvS2v0+XqMl/FTH9hGai/uRJE8MvWeM8yVdYft9zX3rc6Hth4eddFC7JO1K8r9nRFs0\nF/xSdZGk95LMJvlG0hOSzht404KGCvxlSafaPtn2kZp7oeLJgbYsyLY19z3iTJK7ht4zTpLbk6xJ\nsk5zj+vzSZbkVSbJJ5J22j5t/lMbJL054KRxPpB0ju0V818XG7SEXxSU5p4i/eiSfGv7eklbNfdK\n5INJdgyxpYPzJV0r6XXbr81/7s9Jnh5wUyU3SHpk/j/070q6buA9B5Vkm+0tkrZr7qcrr2qJv6uN\nd7IBhfFONqAwAgcKI3CgMAIHCiNwoDACBwojcKAwAgcK+y+Qlolk5IOqTAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xb0de518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import os\n",
    "import tempfile\n",
    "import subprocess\n",
    "import tempfile\n",
    "import shutil\n",
    "\n",
    "import numpy as np\n",
    "import nibabel as nib\n",
    "os.environ[\"FSLOUTPUTTYPE\"] = \"NIFTI_GZ\"\n",
    "\n",
    "def binarise(data):\n",
    "    # Remember the directory where we started\n",
    "    cwd_orig = os.getcwd()\n",
    "    try:\n",
    "        # Create a temporary directory\n",
    "        tempdir = tempfile.mkdtemp(\"fsl\")\n",
    "        \n",
    "        # Save input data in temp directory\n",
    "        os.chdir(tempdir)\n",
    "        tmpin = nib.Nifti1Image(data, np.identity(4))\n",
    "        tmpin.to_filename(\"in.nii.gz\")\n",
    "        \n",
    "        # Run a command from $FSLDIR\n",
    "        fslmaths = os.path.join(os.environ[\"FSLDIR\"], \"bin\", \"fslmaths\")\n",
    "        \n",
    "        # We could use os.system here if we don't care about returning the stdout/stderr\n",
    "        p = subprocess.Popen([fslmaths, \"in.nii.gz\", \"-thr\", \"0.5\", \"-bin\", \"out.nii.gz\"], \n",
    "                             stdout=subprocess.PIPE, stderr=subprocess.STDOUT)\n",
    "        cmd_stdout = \"\"\n",
    "        while 1:\n",
    "            retcode = p.poll()\n",
    "            cmd_stdout += p.stdout.readline()\n",
    "            if retcode is not None: break\n",
    "        if retcode != 0:\n",
    "            raise RuntimeError(\"Error: %s\" % cmd_stdout)\n",
    "        \n",
    "        # Load the output file and return it with the command standard output\n",
    "        out_nii = nib.load(\"out.nii.gz\")\n",
    "        return out_nii.get_data(), cmd_stdout\n",
    "    finally:\n",
    "        # Change back to our starting directory\n",
    "        os.chdir(cwd_orig)\n",
    "\n",
    "data = np.random.rand(10, 10, 10, 10)\n",
    "plt.imshow(data[:,:,5,5])\n",
    "plt.show()\n",
    "\n",
    "output, stdout = binarise(data)\n",
    "plt.imshow(output[:,:,5,5])\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary\n",
    "\n",
    "## What we've done\n",
    "\n",
    " - It's not that hard to call existing C++ code from Python\n",
    " - Need to be a bit careful with Numpy arrays\n",
    " - Cython is probably the easiest method\n",
    " - Data copying can be minimised by passing data as C arrays (`float *` etc)\n",
    " - `ctypes` may be a good alternative if you have binary compatibility issues\n",
    " - Can always wrap a command line tool as a way of getting started!\n",
    " \n",
    " \n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}