#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION

#include <stdlib.h>
#include <Python.h>
#include <numpy/arrayobject.h>
#include "cshape.h"

static char module_docstring[] = "This module links to C-compiled code for efficient calculation of the surface area "
                                 "in the pyRadiomics package. It provides fast calculation using a marching cubes "
                                 "algortihm, accessed via ""calculate_surfacearea"". Arguments for this function"
                                 "are positional and consist of two numpy arrays, mask and pixelspacing, which must "
                                 "be supplied in this order. Pixelspacing is a 3 element vector containing the pixel"
                                 "spacing in z, y and x dimension, respectively. All non-zero elements in mask are "
                                 "considered to be a part of the segmentation and are included in the algorithm.";
static char coefficients_docstring[] = "Arguments: Mask, PixelSpacing. Uses a marching cubes algorithm to calculate an "
                                       "approximation to the total surface area, volume and maximum diameters. "
                                       "The isovalue is considered to be situated midway between a voxel that is part "
                                       "of the segmentation and a voxel that is not.";
static char coefficients2D_docstring[] = "Arguments: Mask, PixelSpacing. Uses an adapted 2D marching cubes algorithm "
                                         "to calculate an approximation to the total perimeter, surface and maximum "
                                         "diameter. The isovalue is considered to be situated midway between a pixel "
                                         "that is part of the segmentation and a pixel that is not.";

static PyObject *cshape_calculate_coefficients(PyObject *self, PyObject *args);
static PyObject *cshape_calculate_coefficients2D(PyObject *self, PyObject *args);

int check_arrays(PyArrayObject *mask_arr, PyArrayObject *spacing_arr, int *size, int *strides, int dimension);

static PyMethodDef module_methods[] = {
  //{"calculate_", cmatrices_, METH_VARARGS, _docstring},
  { "calculate_coefficients", cshape_calculate_coefficients, METH_VARARGS, coefficients_docstring },
  { "calculate_coefficients2D", cshape_calculate_coefficients2D, METH_VARARGS, coefficients2D_docstring },
  { NULL, NULL, 0, NULL }
};

#if PY_MAJOR_VERSION >= 3

static struct PyModuleDef moduledef = {
  PyModuleDef_HEAD_INIT,
  "_cshape",           /* m_name */
  module_docstring,    /* m_doc */
  -1,                  /* m_size */
  module_methods,      /* m_methods */
  NULL,                /* m_reload */
  NULL,                /* m_traverse */
  NULL,                /* m_clear */
  NULL,                /* m_free */
};

#endif

static PyObject *
moduleinit(void)
{
    PyObject *m;

#if PY_MAJOR_VERSION >= 3
    m = PyModule_Create(&moduledef);
#else
    m = Py_InitModule3("_cshape",
                       module_methods, module_docstring);
#endif

  if (m == NULL)
      return NULL;

  return m;
}

#if PY_MAJOR_VERSION < 3
  PyMODINIT_FUNC
  init_cshape(void)
  {
    // Initialize numpy functionality
    import_array();

    moduleinit();
  }
#else
  PyMODINIT_FUNC
  PyInit__cshape(void)
  {
    // Initialize numpy functionality
    import_array();

    return moduleinit();
  }
#endif

static PyObject *cshape_calculate_coefficients(PyObject *self, PyObject *args)
{
  PyObject *mask_obj, *spacing_obj;
  PyArrayObject *mask_arr, *spacing_arr;
  int size[3];
  int strides[3];
  char *mask;
  double *spacing;
  double SA, Volume;
  double diameters[4];
  PyObject *diameter_obj;

  // Parse the input tuple
  if (!PyArg_ParseTuple(args, "OO", &mask_obj, &spacing_obj))
    return NULL;

  // Interpret the input as numpy arrays
  mask_arr = (PyArrayObject *)PyArray_FROM_OTF(mask_obj, NPY_BYTE, NPY_ARRAY_FORCECAST | NPY_ARRAY_IN_ARRAY);
  spacing_arr = (PyArrayObject *)PyArray_FROM_OTF(spacing_obj, NPY_DOUBLE, NPY_ARRAY_FORCECAST | NPY_ARRAY_IN_ARRAY);

  if (check_arrays(mask_arr, spacing_arr, size, strides, 3) > 0) return NULL;

  // Get arrays in Ctype
  mask = (char *)PyArray_DATA(mask_arr);
  spacing = (double *)PyArray_DATA(spacing_arr);

  //Calculate Surface Area and volume
  if (calculate_coefficients(mask, size, strides, spacing, &SA, &Volume, diameters))
  {
    // An error has occurred
    Py_XDECREF(mask_arr);
    Py_XDECREF(spacing_arr);
    PyErr_SetString(PyExc_RuntimeError, "Calculation of Shape coefficients failed.");
    return NULL;
  }

  // Clean up
  Py_XDECREF(mask_arr);
  Py_XDECREF(spacing_arr);

  diameter_obj = Py_BuildValue("ffff", diameters[0], diameters[1], diameters[2], diameters[3]);
  return Py_BuildValue("ffN", SA, Volume, diameter_obj);
}

static PyObject *cshape_calculate_coefficients2D(PyObject *self, PyObject *args)
{
  PyObject *mask_obj, *spacing_obj;
  PyArrayObject *mask_arr, *spacing_arr;
  int size[2];
  int strides[2];
  char *mask;
  double *spacing;
  double perimeter, surface, diameter;

  // Parse the input tuple
  if (!PyArg_ParseTuple(args, "OO", &mask_obj, &spacing_obj))
    return NULL;

  // Interpret the input as numpy arrays
  mask_arr = (PyArrayObject *)PyArray_FROM_OTF(mask_obj, NPY_BYTE, NPY_ARRAY_FORCECAST | NPY_ARRAY_IN_ARRAY);
  spacing_arr = (PyArrayObject *)PyArray_FROM_OTF(spacing_obj, NPY_DOUBLE, NPY_ARRAY_FORCECAST | NPY_ARRAY_IN_ARRAY);

  if (check_arrays(mask_arr, spacing_arr, size, strides, 2) > 0) return NULL;

  // Get arrays in Ctype
  mask = (char *)PyArray_DATA(mask_arr);
  spacing = (double *)PyArray_DATA(spacing_arr);

  //Calculate Surface Area and volume
  if (calculate_coefficients2D(mask, size, strides, spacing, &perimeter, &surface, &diameter))
  {
    // An error has occurred
    Py_XDECREF(mask_arr);
    Py_XDECREF(spacing_arr);
    PyErr_SetString(PyExc_RuntimeError, "Calculation of Shape coefficients failed.");
    return NULL;
  }

  // Clean up
  Py_XDECREF(mask_arr);
  Py_XDECREF(spacing_arr);

  return Py_BuildValue("fff", perimeter, surface, diameter);
}

int check_arrays(PyArrayObject *mask_arr, PyArrayObject *spacing_arr, int *size, int *strides, int dimension)
{
  int i;

  if (mask_arr == NULL || spacing_arr == NULL)
  {
    Py_XDECREF(mask_arr);
    Py_XDECREF(spacing_arr);
    PyErr_SetString(PyExc_RuntimeError, "Error parsing array arguments.");
    return 1;
  }

  if (PyArray_NDIM(mask_arr) != dimension || PyArray_NDIM(spacing_arr) != 1)
  {
    Py_XDECREF(mask_arr);
    Py_XDECREF(spacing_arr);
    PyErr_Format(PyExc_ValueError, "Expected a %iD array for mask, 1D for spacing.", dimension);
    return 2;
  }

  if ( !PyArray_IS_C_CONTIGUOUS(mask_arr) || !PyArray_IS_C_CONTIGUOUS(spacing_arr))
  {
    Py_XDECREF(mask_arr);
    Py_XDECREF(spacing_arr);
    PyErr_SetString(PyExc_ValueError, "Expecting input arrays to be C-contiguous.");
    return 3;
  }

  if (PyArray_DIM(spacing_arr, 0) != PyArray_NDIM(mask_arr))
  {
    Py_XDECREF(mask_arr);
    Py_XDECREF(spacing_arr);
    PyErr_SetString(PyExc_ValueError, "Expecting spacing array to have shape (3,).");
    return 4;
  }

  // Get sizes and strides of the arrays
  for (i = 0; i < dimension; i++)
  {
    size[i] = (int)PyArray_DIM(mask_arr, i);
    strides[i] = (int)(PyArray_STRIDE(mask_arr, i) / PyArray_ITEMSIZE(mask_arr));
  }

  return 0;
}