From 244da2d0d7cc97f6ed616820e30501a9f3a49d8f Mon Sep 17 00:00:00 2001 From: Joost van Griethuysen Date: Mon, 21 Nov 2016 23:19:12 -0500 Subject: [PATCH] STYL: Add docstrings to _cmatrices and _cshape Add additional comments where necessary, explaining the source code. Fix a datatype bug in the python calculation of surface area (must be in a number-like datatype, otherwise the interpolate function will not work properly). --- radiomics/shape.py | 1 + radiomics/src/_cmatrices.c | 20 ++-- radiomics/src/_cshape.c | 15 ++- radiomics/src/cmatrices.c | 1 - radiomics/src/cshape.c | 181 ++++++++++++++++++++----------------- 5 files changed, 122 insertions(+), 96 deletions(-) diff --git a/radiomics/shape.py b/radiomics/shape.py index d4865ad8..00074459 100644 --- a/radiomics/shape.py +++ b/radiomics/shape.py @@ -44,6 +44,7 @@ def __init__(self, inputImage, inputMask, **kwargs): self.SurfaceArea = self._calculateSurfaceArea() def _calculateSurfaceArea(self): + self.maskArray = self.maskArray.astype('int') # needed for the interpolate function to work correctly # define relative locations of the 8 voxels of a sampling cube gridAngles = numpy.array([(0, 0, 0), (0, 0, 1), (0, 1, 1), (0, 1, 0), (1, 0, 0), (1, 0, 1), (1, 1, 1), (1, 1, 0)]) diff --git a/radiomics/src/_cmatrices.c b/radiomics/src/_cmatrices.c index 250da3b5..f766791c 100644 --- a/radiomics/src/_cmatrices.c +++ b/radiomics/src/_cmatrices.c @@ -2,12 +2,18 @@ #include #include "cmatrices.h" -static char module_docstring[] = ""; -static char glcm_docstring[] = ""; -static char gldm_docstring[] = ""; -static char ngtdm_docstring[] = ""; -static char glszm_docstring[] = ""; -static char glrlm_docstring[] = ""; +static char module_docstring[] = ("This module links to C-compiled code for efficient calculation of various matrices " + "in the pyRadiomics package. It provides fast calculation for GLCM, GLDM, NGTDM, " + "GLRLM and GLSZM. All functions are given names as follows: ""calculate_"", " + "where is the name of the matix, in lowercase. Arguments for these functions " + "are positional and start with 3 numpy arrays (image, mask and angles) and 1 integer " + "(Ng, number of gray levels). Optionally extra arguments may be required, see function " + "docstrings for detailed information."); +static char glcm_docstring[] = "Arguments: Image, Mask, Angles, Ng."; +static char gldm_docstring[] = "Arguments: Image, Mask, Angles, Ng, Alpha."; +static char ngtdm_docstring[] = "Arguments: Image, Mask, Angles, Ng."; +static char glszm_docstring[] = "Arguments: Image, Mask, Angles, Ng, Ns, matrix is cropped to maximum size encountered."; +static char glrlm_docstring[] = "Arguments: Image, Mask, Angles, Ng, Nr."; static PyObject *cmatrices_calculate_glcm(PyObject *self, PyObject *args); static PyObject *cmatrices_calculate_gldm(PyObject *self, PyObject *args); @@ -21,7 +27,7 @@ static PyMethodDef module_methods[] = { { "calculate_gldm", cmatrices_calculate_gldm, METH_VARARGS, gldm_docstring }, { "calculate_ngtdm", cmatrices_calculate_ngtdm, METH_VARARGS, ngtdm_docstring }, { "calculate_glszm", cmatrices_calculate_glszm, METH_VARARGS, glszm_docstring }, - { "calculate_glrlm", cmatrices_calculate_glrlm, METH_VARARGS, glszm_docstring }, + { "calculate_glrlm", cmatrices_calculate_glrlm, METH_VARARGS, glrlm_docstring }, { NULL, NULL, 0, NULL } }; diff --git a/radiomics/src/_cshape.c b/radiomics/src/_cshape.c index ded40db8..d04e854b 100644 --- a/radiomics/src/_cshape.c +++ b/radiomics/src/_cshape.c @@ -2,8 +2,16 @@ #include #include "cshape.h" -static char module_docstring[] = ""; -static char surface_docstring[] = ""; +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 surface_docstring[] = "Arguments: Mask, PixelSpacing, uses a marching cubes algorithm to calculate an " + "approximation to the total surface area. The isovalue is considered to be situated " + "midway between a voxel that is part of the segmentation and a voxel that is not."; static PyObject *cshape_calculate_surfacearea(PyObject *self, PyObject *args); @@ -40,6 +48,7 @@ static PyObject *cshape_calculate_surfacearea(PyObject *self, PyObject *args) mask_arr = (PyArrayObject *)PyArray_FROM_OTF(mask_obj, NPY_BYTE, NPY_FORCECAST | NPY_UPDATEIFCOPY | NPY_IN_ARRAY); spacing_arr = (PyArrayObject *)PyArray_FROM_OTF(spacing_obj, NPY_DOUBLE, NPY_FORCECAST | NPY_UPDATEIFCOPY | NPY_IN_ARRAY); + // Check if there were any errors during casting to array objects if (mask_arr == NULL || spacing_arr == NULL) { Py_XDECREF(mask_arr); @@ -72,7 +81,7 @@ static PyObject *cshape_calculate_surfacearea(PyObject *self, PyObject *args) Py_DECREF(mask_arr); Py_DECREF(spacing_arr); - if (SA < 0) + if (SA < 0) // if SA < 0, an error has occurred { PyErr_SetString(PyExc_RuntimeError, "Calculation of GLCM Failed."); return NULL; diff --git a/radiomics/src/cmatrices.c b/radiomics/src/cmatrices.c index 3f3cc7fb..b54f4364 100644 --- a/radiomics/src/cmatrices.c +++ b/radiomics/src/cmatrices.c @@ -1,4 +1,3 @@ -#include #include "cmatrices.h" int calculate_glcm(int *image, char *mask, int Sx, int Sy, int Sz, int *angles, int Na, double *glcm, int Ng) diff --git a/radiomics/src/cshape.c b/radiomics/src/cshape.c index 1959d05b..6e602e91 100644 --- a/radiomics/src/cshape.c +++ b/radiomics/src/cshape.c @@ -1,6 +1,101 @@ #include "cmatrices.h" #include +// Declare the look-up tables, these are filled at the bottom of this code file. +static const int gridAngles[8][3]; +static const int edgeTable[128]; +static const int triTable[128][16]; + +double calculate_surfacearea(char *mask, int *size, double *spacing) +{ + int iz, iy, ix, i, j, t, d; // iterator indices + unsigned char cube_idx; // cube identifier, 8 bits signifying which corners of the cube belong to the segmentation + int ang; // Angle index (8 'angles', one pointing to each corner of the marching cube + double vertList[12][3]; // array to hold location of isosurface vectors (points) on the marching cube edges. + double sum; + double surfaceArea = 0; // Total surface area + double a[3], b[3], c[3]; // 2 points of the triangle, relative to the third, and the cross product vector + + // Iterate over all voxels, do not include last voxels in the three dimensions, as the cube includes voxels at pos +1 + for (iz = 0; iz < (size[0] - 1); iz++) + { + for (iy = 0; iy < (size[1] - 1); iy++) + { + for (ix = 0; ix < (size[2] - 1); ix++) + { + // Get current cube_idx by analyzing each point of the current cube + cube_idx = 0; + i = ix + iy * size[2] + iz * size[2] * size[1]; // voxel at cube corner (0, 0, 0) + for (ang = 0; ang < 8; ang++) + { + j = i + gridAngles[ang][0] * size[2] * size[1] + gridAngles[ang][1] * size[2] + gridAngles[ang][2]; + if (mask[j]) cube_idx |= (1 << ang); + } + // Isosurface is symmetrical around the midpoint, flip the number if > 128 + // This enables look-up tables to be 1/2 the size. + if (cube_idx & 0x80) cube_idx ^= 0xff; + + // Exlcude cubes entirely outside or inside the segmentation. Also exclude potential invalid values (> 128). + if (cube_idx > 0 && cube_idx < 128) + { + if (edgeTable[cube_idx] & 1) interpolate(vertList[0], 1, 0, spacing); + if (edgeTable[cube_idx] & 2) interpolate(vertList[1], 2, 1, spacing); + if (edgeTable[cube_idx] & 4) interpolate(vertList[2], 2, 3, spacing); + if (edgeTable[cube_idx] & 8) interpolate(vertList[3], 3, 0, spacing); + if (edgeTable[cube_idx] & 16) interpolate(vertList[4], 5, 4, spacing); + if (edgeTable[cube_idx] & 32) interpolate(vertList[5], 6, 5, spacing); + if (edgeTable[cube_idx] & 64) interpolate(vertList[6], 6, 7, spacing); + if (edgeTable[cube_idx] & 128) interpolate(vertList[7], 7, 4, spacing); + if (edgeTable[cube_idx] & 256) interpolate(vertList[8], 4, 0, spacing); + if (edgeTable[cube_idx] & 512) interpolate(vertList[9], 5, 1, spacing); + if (edgeTable[cube_idx] & 1024) interpolate(vertList[10], 6, 2, spacing); + if (edgeTable[cube_idx] & 2048) interpolate(vertList[11], 7, 3, spacing); + + t = 0; + while (triTable[cube_idx][t*3] >= 0) // Exit loop when no more triangles are present (element at index = -1) + { + for (d = 0; d < 3; d++) + { + a[d] = vertList[triTable[cube_idx][t*3 + 1]][d] - vertList[triTable[cube_idx][t*3]][d]; + b[d] = vertList[triTable[cube_idx][t*3 + 2]][d] - vertList[triTable[cube_idx][t*3]][d]; + } + c[0] = (a[1] * b[2]) - (b[1] * a[2]); + c[1] = (a[2] * b[0]) - (b[2] * a[0]); + c[2] = (a[0] * b[1]) - (b[0] * a[1]); + + // Get the square + c[0] = c[0] * c[0]; + c[1] = c[1] * c[1]; + c[2] = c[2] * c[2]; + + // Compute the surface, which is equal to 1/2 magnitude of the cross product, where + // The magnitude is obtained by calculating the euclidean distance between (0, 0, 0) + // and the location of c + sum = c[0] + c[1] + c[2]; + sum = sqrt(sum); + sum = 0.5 * sum; + surfaceArea += sum; + t++; + } + } + } + } + } + return surfaceArea; +} + +int interpolate(double *vertEntry, int a1, int a2, double *spacing) +{ + int d; + + for (d = 0; d < 3; d++) + { + if (gridAngles[a1][d] == 1 && gridAngles[a2][d] == 1) vertEntry[d] = spacing[d]; + else if (gridAngles[a1][d] != gridAngles[a2][d] ) vertEntry[d] = 0.5 * spacing[d]; + else vertEntry[d] = 0; + } +} + static const int gridAngles[8][3] = { { 0, 0, 0 }, { 0, 0, 1 }, { 0, 1, 1 }, {0, 1, 0}, { 1, 0, 0 }, {1, 0, 1 }, { 1, 1, 1 }, { 1, 1, 0 }}; static const int edgeTable[128] = { 0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, @@ -141,88 +236,4 @@ static const int triTable[128][16] = { { 0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }, { 7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1 }, { 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } -}; - -double calculate_surfacearea(char *mask, int *size, double *spacing) -{ - int iz, iy, ix, i, j, t, d; - unsigned char cube_idx; - int ang; - int Ni = size[0] * size[1] * size[2]; - double sum; - double SA = 0; - double vertList[12][3]; - double a[3], b[3], c[3]; - - for (iz = 0; iz < (size[0] - 1); iz++) - { - for (iy = 0; iy < (size[1] - 1); iy++) - { - for (ix = 0; ix < (size[2] - 1); ix++) - { - cube_idx = 0; - i = ix + iy * size[2] + iz * size[2] * size[1]; - for (ang = 0; ang < 8; ang++) - { - j = i + gridAngles[ang][0] * size[2] * size[1] + gridAngles[ang][1] * size[2] + gridAngles[ang][2]; - if (mask[j]) cube_idx |= (1 << ang); - } - if (cube_idx & 0x80) cube_idx ^= 0xff; - if (cube_idx > 0 && cube_idx < 128) - { - if (edgeTable[cube_idx] & 1) interpolate(vertList[0], 1, 0, spacing); - if (edgeTable[cube_idx] & 2) interpolate(vertList[1], 2, 1, spacing); - if (edgeTable[cube_idx] & 4) interpolate(vertList[2], 2, 3, spacing); - if (edgeTable[cube_idx] & 8) interpolate(vertList[3], 3, 0, spacing); - if (edgeTable[cube_idx] & 16) interpolate(vertList[4], 5, 4, spacing); - if (edgeTable[cube_idx] & 32) interpolate(vertList[5], 6, 5, spacing); - if (edgeTable[cube_idx] & 64) interpolate(vertList[6], 6, 7, spacing); - if (edgeTable[cube_idx] & 128) interpolate(vertList[7], 7, 4, spacing); - if (edgeTable[cube_idx] & 256) interpolate(vertList[8], 4, 0, spacing); - if (edgeTable[cube_idx] & 512) interpolate(vertList[9], 5, 1, spacing); - if (edgeTable[cube_idx] & 1024) interpolate(vertList[10], 6, 2, spacing); - if (edgeTable[cube_idx] & 2048) interpolate(vertList[11], 7, 3, spacing); - - for (t = 0; t < 5; t++) - { - if (triTable[cube_idx][t*3] >= 0) - { - for (d = 0; d < 3; d++) - { - a[d] = vertList[triTable[cube_idx][t*3 + 1]][d] - vertList[triTable[cube_idx][t*3]][d]; - b[d] = vertList[triTable[cube_idx][t*3 + 2]][d] - vertList[triTable[cube_idx][t*3]][d]; - } - c[0] = (a[1] * b[2]) - (b[1] * a[2]); - c[1] = (a[2] * b[0]) - (b[2] * a[0]); - c[2] = (a[0] * b[1]) - (b[0] * a[1]); - - // Get the square - c[0] = c[0] * c[0]; - c[1] = c[1] * c[1]; - c[2] = c[2] * c[2]; - - sum = c[0] + c[1] + c[2]; - sum = sqrt(sum); - sum = 0.5 * sum; - SA += sum; - } - } - } - } - } - } - return SA; -} - -int interpolate(double *vertEntry, int a1, int a2, double *spacing) -{ - int d; - - for (d = 0; d < 3; d++) - { - if (gridAngles[a1][d] == 1 && gridAngles[a2][d] == 1) vertEntry[d] = spacing[d]; - else if (gridAngles[a1][d] == 0 && gridAngles[a2][d] == 1) vertEntry[d] = 0.5 * spacing[d]; - else if (gridAngles[a1][d] == 1 && gridAngles[a2][d] == 0) vertEntry[d] = 0.5 * spacing[d]; - else vertEntry[d] = 0; - } -} \ No newline at end of file +}; \ No newline at end of file