ENH: Import the ITKBridgeNumPy module

This is imported from https://github.com/InsightSoftwareConsortium/ITKBridgeNumPy.git
@e74bdafc2824c7f63fdc7c0e36ae843db0c92393

This Remote Module is mature and is always used when building the Python
wrapping. Import it into the repository to make it easier for packagers,
so they do not have to package it separately or worry about pulling it
down at CMake configure time.

Change-Id: Ic338eebb6502c8876fb8c0707e1a2e0623223043

CMakeLists.txt

@@ -211,10 +211,6 @@ if(ITK_WRAP_PYTHON OR ITK_WRAP_JAVA OR ITK_WRAP_RUBY OR ITK_WRAP_PERL OR ITK_WRA
 else()
   set(ITK_WRAPPING OFF CACHE INTERNAL "Build external languages support" FORCE)
 endif()
-# Enable BridgeNumPy by default when Python wrapping is enabled
-if(ITK_WRAP_PYTHON AND NOT DEFINED Module_BridgeNumPy)
-  set(Module_BridgeNumPy ON CACHE BOOL "Enable the BridgeNumPy module.")
-endif()
 
 include(ITKSetStandardCompilerFlags)
 #---------------------------------------------------------------

Modules/Bridge/NumPy/CMakeLists.txt

@@ -0,0 +1,12 @@
+cmake_minimum_required(VERSION 3.9.5)
+project(BridgeNumPy)
+
+set(BridgeNumPy_SYSTEM_INCLUDE_DIRS "${PYTHON_INCLUDE_DIR}")
+
+if(NOT ITK_SOURCE_DIR)
+  find_package(ITK REQUIRED)
+  list(APPEND CMAKE_MODULE_PATH ${ITK_CMAKE_DIR})
+  include(ITKModuleExternal)
+else()
+  itk_module_impl()
+endif()

Modules/Bridge/NumPy/README.rst

@@ -0,0 +1,82 @@
+ITKBridgeNumPy
+==============
+
+This is a port from the original WrapITK PyBuffer to an ITKv4 module.
+
+Differences from the original PyBuffer:
+
+- Support for VectorImage's
+- Option to not swap the axes (only for GetArrayViewFromImage)
+- Tests
+- Based on the `Python Buffer Protocol <https://docs.python.org/3/c-api/buffer.html>`_ -- only requires NumPy at run time, not build time.
+
+This module is available in the ITK source tree as a Remote
+module.  To enable it, set::
+
+  ITK_WRAP_PYTHON:BOOL=ON
+  Module_ITKBridgeNumPy:BOOL=ON
+
+in ITK's CMake build configuration. In ITK 4.12 and later, this module is
+enabled by default when `ITK_WRAP_PYTHON` is enabled.
+
+To get a view of an ITK image in a NumPy array::
+
+  import itk
+
+  PixelType = itk.ctype('float')
+  Dimension = 3
+  ImageType = itk.Image[PixelType, Dimension]
+
+  image = ImageType.New()
+  size = itk.Size[Dimension]()
+  size.Fill(100)
+  region = itk.ImageRegion[Dimension](size)
+  image.SetRegions(region)
+  image.Allocate()
+
+  arr = itk.PyBuffer[ImageType].GetArrayViewFromImage(image)
+
+To get a view of a NumPy array in an ITK image::
+
+  import numpy as np
+  import itk
+
+  PixelType = itk.ctype('float')
+  Dimension = 3
+  ImageType = itk.Image[PixelType, Dimension]
+
+  arr = np.zeros((100, 100, 100), np.float32)
+  image = itk.PyBuffer[ImageType].GetImageViewFromArray(arr)
+
+It is also possible to get views of VNL matrices and arrays from NumPy arrays and
+back::
+
+  import numpy as np
+  import itk
+
+  ElementType = itk.ctype('float')
+  vector = itk.vnl_vector[ElementType]()
+  vector.set_size(8)
+  arr = itk.PyVnl[ElementType].GetArrayViewFromVnlVector(vector)
+
+  matrix = itk.vnl_matrix[ElementType]()
+  matrix.set_size(3, 4)
+  arr = itk.PyVnl[ElementType].GetArrayViewFromVnlMatrix(matrix)
+
+  arr = np.zeros((100,), np.float32)
+  vector = itk.PyVnl[ElementType].GetVnlVectorViewFromArray(arr)
+
+  arr = np.zeros((100, 100), np.float32)
+  matrix = itk.PyVnl[ElementType].GetVnlMatrixViewFromArray(arr)
+
+.. warning::
+
+  The conversions create `NumPy Views
+  <https://scipy-cookbook.readthedocs.io/items/ViewsVsCopies.html>`_, i.e. it
+  presents the ITK image pixel buffer in the NumPy array, and the buffer is
+  shared. This means that no copies are made (which increases speed and
+  reduces memory consumption). It also means that any changes in the NumPy
+  array change the ITK image content. Additionally, a reference to an ITK
+  image object must be available to use its NumPy array view. Using an array
+  view after its source image has been deleted can results in corrupt values
+  or a segfault.

Modules/Bridge/NumPy/include/itkPyBuffer.h

@@ -0,0 +1,94 @@
+/*=========================================================================
+ *
+ *  Copyright Insight Software Consortium
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+
+#ifndef itkPyBuffer_h
+#define itkPyBuffer_h
+
+#include "itkObject.h"
+#include "itkObjectFactory.h"
+#include "itkImportImageFilter.h"
+#include "itkDefaultConvertPixelTraits.h"
+
+
+// The python header defines _POSIX_C_SOURCE without a preceding #undef
+#undef _POSIX_C_SOURCE
+#undef _XOPEN_SOURCE
+// For Python 2.7 hypot bug, see https://bugs.python.org/issue11566
+#include "PatchedPython27pyconfig.h"
+#include <Python.h>
+
+namespace itk
+{
+
+/** \class PyBuffer
+ *
+ *  \brief Helper class to get ITK image views into python arrays and back.
+ *
+ *  This class will receive a C buffer and create the equivalent python
+ *  array view. This permits passing image buffers into python arrays from
+ *  the NumPy python package.
+ *
+ *  \ingroup ITKBridgeNumPy
+ */
+template <typename TImage>
+class PyBuffer
+{
+public:
+  /** Standard "Self" typedef. */
+  typedef PyBuffer         Self;
+
+  /** Type of the image from which the buffer will be converted */
+  typedef TImage                                                       ImageType;
+  typedef typename ImageType::PixelType                                PixelType;
+  typedef typename ImageType::SizeType                                 SizeType;
+  typedef typename ImageType::IndexType                                IndexType;
+  typedef typename ImageType::RegionType                               RegionType;
+  typedef typename ImageType::PointType                                PointType;
+  typedef typename ImageType::SpacingType                              SpacingType;
+  typedef typename ImageType::Pointer                                  ImagePointer;
+  typedef typename DefaultConvertPixelTraits<PixelType>::ComponentType ComponentType;
+
+   /** Image dimension. */
+  itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension);
+
+  typedef typename Image< ComponentType, ImageDimension >::Pointer  OutputImagePointer;
+
+  /**
+   * Get an Array with the content of the image buffer
+   */
+  static PyObject * _GetArrayViewFromImage( ImageType * image);
+
+  /**
+   * Get an ITK image from a Python array
+   */
+  static const OutputImagePointer _GetImageViewFromArray( PyObject *arr, PyObject *shape, PyObject *numOfComponent);
+
+protected:
+
+private:
+  PyBuffer(const Self&);       // Purposely not implemented.
+  void operator=(const Self&); // Purposely not implemented.
+};
+
+} // namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkPyBuffer.hxx"
+#endif
+
+#endif // _itkPyBuffer_h

Modules/Bridge/NumPy/include/itkPyBuffer.hxx

@@ -0,0 +1,182 @@
+/*=========================================================================
+ *
+ *  Copyright Insight Software Consortium
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkPyBuffer_hxx
+#define itkPyBuffer_hxx
+
+#include "itkPyBuffer.h"
+
+namespace itk
+{
+
+template<class TImage>
+PyObject *
+PyBuffer<TImage>
+::_GetArrayViewFromImage( ImageType * image)
+{
+  PyObject *                  memoryView    = NULL;
+  Py_buffer                   pyBuffer;
+  memset(&pyBuffer, 0, sizeof(Py_buffer));
+
+  Py_ssize_t                  len           = 1;
+  size_t                      pixelSize     = sizeof(ComponentType);
+  int                         res           = 0;
+
+  if( !image )
+    {
+    throw std::runtime_error("Input image is null");
+    }
+
+  image->Update();
+
+  ComponentType *buffer =  const_cast < ComponentType *> (reinterpret_cast< const ComponentType* > ( image->GetBufferPointer() ) );
+
+  void * itkImageBuffer = (void *)( buffer );
+
+  // Computing the length of data
+  const int numberOfComponents = image->GetNumberOfComponentsPerPixel();
+  SizeType size = image->GetBufferedRegion().GetSize();
+
+  for( unsigned int dim = 0; dim < ImageDimension; ++dim )
+    {
+    len *= size[dim];
+    }
+
+  len *= numberOfComponents;
+  len *= pixelSize;
+
+  res = PyBuffer_FillInfo(&pyBuffer, NULL, (void*)itkImageBuffer, len, 0, PyBUF_CONTIG);
+  memoryView = PyMemoryView_FromBuffer(&pyBuffer);
+
+  PyBuffer_Release(&pyBuffer);
+
+  return memoryView;
+}
+
+template<class TImage>
+const typename PyBuffer<TImage>::OutputImagePointer
+PyBuffer<TImage>
+::_GetImageViewFromArray( PyObject *arr, PyObject *shape, PyObject *numOfComponent)
+{
+  PyObject *                  shapeseq      = NULL;
+  PyObject *                  item          = NULL;
+
+  Py_ssize_t                  bufferLength;
+  Py_buffer                   pyBuffer;
+  memset(&pyBuffer, 0, sizeof(Py_buffer));
+
+  SizeType size;
+  SizeType sizeFortran;
+  SizeValueType numberOfPixels = 1;
+
+  const void *                buffer;
+
+  long                        numberOfComponents= 1;
+  unsigned int                dimension     = 0;
+
+
+  size_t                      pixelSize     = sizeof(ComponentType);
+  size_t                      len           = 1;
+
+  if(PyObject_GetBuffer(arr, &pyBuffer, PyBUF_WRITABLE | PyBUF_ND | PyBUF_ANY_CONTIGUOUS ) == -1)
+    {
+    PyErr_SetString( PyExc_RuntimeError, "Cannot get an instance of NumPy array." );
+    PyBuffer_Release(&pyBuffer);
+    return NULL;
+    }
+  else
+    {
+    bufferLength = pyBuffer.len;
+    buffer = pyBuffer.buf;
+    }
+  PyBuffer_Release(&pyBuffer);
+
+  shapeseq   = PySequence_Fast(shape, "expected sequence");
+  dimension  = PySequence_Size(shape);
+
+  numberOfComponents = PyInt_AsLong(numOfComponent);
+
+  for( unsigned int i = 0; i < dimension; ++i )
+    {
+    item = PySequence_Fast_GET_ITEM(shapeseq,i);
+    size[i] = (SizeValueType)PyInt_AsLong(item);
+    sizeFortran[dimension - 1 - i] = (SizeValueType)PyInt_AsLong(item);
+    numberOfPixels *= size[i];
+    }
+
+  bool isFortranContiguous = false;
+  if( pyBuffer.strides != NULL && pyBuffer.itemsize == pyBuffer.strides[0] )
+    {
+    isFortranContiguous = true;
+    }
+
+  len = numberOfPixels*numberOfComponents*pixelSize;
+  if ( bufferLength != len )
+    {
+    PyErr_SetString( PyExc_RuntimeError, "Size mismatch of image and Buffer." );
+    PyBuffer_Release(&pyBuffer);
+    Py_DECREF(shapeseq);
+    return NULL;
+    }
+
+  IndexType start;
+  start.Fill( 0 );
+
+  RegionType region;
+  region.SetIndex( start );
+  region.SetSize( size );
+  if( isFortranContiguous )
+    {
+    region.SetSize( sizeFortran );
+    }
+  else
+    {
+    region.SetSize( size );
+    }
+
+  PointType origin;
+  origin.Fill( 0.0 );
+
+  SpacingType spacing;
+  spacing.Fill( 1.0 );
+
+  typedef ImportImageFilter< ComponentType, ImageDimension > ImporterType;
+  typename ImporterType::Pointer importer = ImporterType::New();
+  importer->SetRegion( region );
+  importer->SetOrigin( origin );
+  importer->SetSpacing( spacing );
+  const bool importImageFilterWillOwnTheBuffer = false;
+
+  ComponentType * data = (ComponentType *)buffer;
+
+  importer->SetImportPointer( data,
+                              numberOfPixels,
+                              importImageFilterWillOwnTheBuffer );
+
+  importer->Update();
+  OutputImagePointer output = importer->GetOutput();
+  output->DisconnectPipeline();
+
+  Py_DECREF(shapeseq);
+  PyBuffer_Release(&pyBuffer);
+
+  return output;
+}
+
+} // namespace itk
+
+#endif