通过带distutils的numpy阵列的元素未对准

Question

我正在CentOS 7上使用Anaconda Spyder 3.3.6（Python 3.7）

我有一个4D numpy数组，surfaceVectors，其中第四维是3D体积中每个体素处的非规范化向量（大小和方向）。我制作了一个对应的3D体积，大小，其中每个元素是surfaceVectors中对应元素处向量的范数。我做的震级如下

dimensions=filteredVolumes[0].shape
numVoxels = int(len(surfaceVectors.reshape(-1))/3)
temp=surfaceVectors.reshape(numVoxels, 3)
temp2=np.array([np.linalg.norm(temp[item],axis=0) for item in range(0,numVoxels)])
magnitudes = temp2.reshape(dimensions)

然后我加载一个模块surfaceModules，其形式如下

#include <stdio.h>
#include <stdlib.h>
#include "/home/peter/anaconda3/include/python3.7m/Python.h"
#include "numpy/arrayobject.h"

#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION

int  not_doublematrix(PyArrayObject *mat);
double ****pymatrix_to_Carrayptrs(PyArrayObject *arrayin);
double **ptrvector(long n, long l);
double ****pptrvector(long n, long m, long l);
void free_Carrayptrs(double ****v, int m, int n, int l);

PyArrayObject * CsqueezeSurfaces(PyArrayObject **surfaceVectors, PyArrayObject **    dimensions)
{
    printf("This function squeezes the surfaces");

    return *surfaceVectors;
}

/* ==== Allocate a double *vector (vec of pointers) ======================
    Memory is Allocated!  See void free_Carray(double ** )                  */
double **ptrvector(long n, long l)  {
    int i;
    double **v;
    if (!(v=(double **)malloc((size_t) (n*sizeof(double *))))){
        printf("In **ptrvector. Allocation of memory for double array failed.");
        return NULL;  }
    for (i=0; i<n; ++i){
        if (!(v[i]=(double *)malloc((size_t) (l*sizeof(double))))){
            for (--i; i>=0; --i) free(v[i]);
        }
        printf("In **ptrvector. Allocation of memory for double array failed.");
        free(v);
        return NULL;
    }
    return v;
}
double ****pptrvector(long n, long m, long l)  {
    int i, j, k;
    double ****v;
    if (!(v=(double ****)malloc((size_t) (n*sizeof(double ***)))))   {
        printf("In ***pptrvector. Allocation of memory for double array failed.");
        return NULL;
    }
    for (i=0; i<n; ++i){
        if (!(v[i]=(double ***)malloc((size_t) (m*sizeof(double **))))){
            for (--i; i>=0; --i){
            free(v[i]);
            }

            free(v);
            printf("In ***pptrvector. Allocation of memory for double array failed.");
            return NULL;
            }
            for (j=0; j<m; ++j){
                if (!(v[i][j]=(double **)malloc((size_t) (l*sizeof(double *))))){
                    for (--j; j>=0; --j) {
                        for (k=0; k<l; ++k) free(v[i][j][k]);
                        free(v[i][j]);
                    }
                    free(v[i]);
                    for (--i; i>=0; --i){
                        for (j=0; j<m; ++j) {
                            for (k=0; k<l; ++k) free(v[i][j][k]);
                            free(v[i][j]);
                        }
                        free(v[i]);
                    }
                    free(v);
                    printf("In ***pptrvector. Allocation of memory for double array failed.");
                    return NULL;
                }
                for (k=0; k<l; ++k){
                    if (!(v[i][j][k]=(double *)malloc((size_t) (3*sizeof(double))))){
                        for (--k; k>=0; --k) free(v[i][j][k]);
                        free(v[i][j]);
                        for (--j;j>=0; --j){
                            for (k=0; k<l; ++k) free(v[i][j][k]);
                            free(v[i][j]);
                        }
                        free(v[i]);
                        for (--i; i>=0; --i) {
                            for (j=0; j<m; ++j){
                                for (k=0; k<l; ++k) free(v[i][j][k]);
                                free(v[i][j]);
                                }
                            free(v[i]);
                        }
                        free(v);
                        printf("In ***pptrvector. Allocation of memory for double array failed.");
                        return NULL;
                }
            }
        }
    }

    return v;
}

void free_Carrayptrs(double ****v, int m, int n, int l) {
    int i, j, k;
    for (i=0; i<m; ++i){
        for (j=0; j<n; ++j){
            for (k=0; k<l; ++k){
                free(v[i][j][k]);
            }
            free(v[i][j]);
        }
        free(v[i]);
    }
    free((char*) v);
}

/* ==== Create Carray from PyArray ======================
    Assumes PyArray is contiguous in memory.
    Memory is allocated!                                    */
double ****pymatrix_to_Carrayptrs(PyArrayObject *arrayin)  {
    double ****d;
    int h,i,j,k,n,m,l, inc=0;

    n=arrayin->dimensions[0];
    m=arrayin->dimensions[1];
    l=arrayin->dimensions[2];
    d=pptrvector(n,m,l);
    for ( i=0; i<n; i++)  {
        for (j=0; j<m; ++j){
            for (h=0; h<l; ++h){
                for (k=0; k<3; ++k)  d[i][j][h][k]=arrayin->data[inc++];
           }
        }
    }

    return d;
}

/* ==== Check that PyArrayObject is a double (Float) type and a matrix ==============
    return 1 if an error and raise exception */
int  not_doublematrix(PyArrayObject *mat)  {
    if (mat->descr->type_num != NPY_DOUBLE || mat->nd != 4)  {
        PyErr_SetString(PyExc_ValueError,
            "In not_doublematrix: array must be of type Float and 4 dimensional (n x m).");
        return 1;  }
    return 0;
}

static PyObject *  squeezeSurfaces(PyObject* self, PyObject* args)
{
    PyArrayObject *matin, *magmat, *matout;
    double ****cin, ****cout, ***cmag;
    int h,i,j,k,n,m,l,index=0, dims[4];

    // Parse tuples separately since args will differ between C fcns
    if (!PyArg_ParseTuple(args, "O!O!",
        &PyArray_Type, &matin, &PyArray_Type, &magmat))  return NULL;
    if (NULL == matin)  return NULL;

    // Check that object input is 'double' type and a matrix
    //   Not needed if python wrapper function checks before call to this routine
    if (not_doublematrix(matin)) return NULL;

    // Get the dimensions of the input
    m=dims[0]=matin->dimensions[0];
    n=dims[1]=matin->dimensions[1];
    l=dims[2]=matin->dimensions[2];
    dims[3]=matin->dimensions[3];

    // Make a new double matrix of same dims
    matout=(PyArrayObject *) PyArray_FromDims(4,dims,NPY_DOUBLE);

    // Change contiguous arrays into C ** arrays (Memory is Allocated!)
    cin=pymatrix_to_Carrayptrs(matin);
    cout=pymatrix_to_Carrayptrs(matout);

    // Do the calculation.
    for ( i=0; i<m; i++)  {
        for ( j=0; j<n; j++)  {
            for (k=0; k<l; ++k){
                if (magmat->data[index]>1){
                    // printf("f(%d,%d,%d) = %d\n",k,j,i,(int)(magmat->data[index]+0.5));
                    printf("f(%d,%d,%d) = %d, v=(%d,%d,%d)\n",k,j,i,(int)(magmat->data[index]+0.5),
                        (int)(cin[i][j][k][0]+0.5),
                        (int)(cin[i][j][k][1]+0.5),(int)(cin[i][j][k][2]+0.5));
                }
                ++index;
            }
    }  
}

    // Free memory, close file and return
    free_Carrayptrs(cin, m, n, l);
    free_Carrayptrs(cout, m, n, l);

    return PyArray_Return(matout);
}

int CmergeDirectionalVolumes(PyArrayObject *directions[3], int dimensions[3], PyArrayObject ***filteredVolumes,
    PyArrayObject ***strongest, PyArrayObject ***secondStrongest, PyArrayObject     ****mergedMagnitudes, PyArrayObject ****mergedDirections)
{
    printf("This function merges 3D directional volumes");

    return 1;
}

static PyObject* mergeDirectionalVolumes(PyObject* self, PyObject* args)
{
    int dimensions[3];
    PyArrayObject *directions[3], ***filteredVolumes=NULL, ***strongest=NULL,     ***secondStrongest=NULL;
    PyArrayObject   ****mergedMagnitudes=NULL, ****mergedDirections=NULL;

    if (!PyArg_ParseTuple(args, "i", directions, dimensions, filteredVolumes, strongest, secondStrongest,
            mergedMagnitudes, mergedDirections))
        return NULL;

        return Py_BuildValue("i", CmergeDirectionalVolumes(directions, dimensions, filteredVolumes, strongest,
            secondStrongest, mergedMagnitudes, mergedDirections));
}

static PyObject* version(PyObject* self){
    return Py_BuildValue("s", "Version 1.0");
}

static PyMethodDef surfaceMethods[] = {
    {"mergeDirectionalVolumes", mergeDirectionalVolumes, METH_VARARGS, "Merges directional volumes"},
    {"squeezeSurfaces", squeezeSurfaces, METH_VARARGS, "Squeezes surfaces along normal vectors"},
    {"version", (PyCFunction)version, METH_NOARGS, "Returns the version"},
    {NULL, NULL, 0, NULL}
};

static struct PyModuleDef surfaceModules = {
    PyModuleDef_HEAD_INIT,
    "surfaceModules",
    "Surface function Module",
    -1,
    surfaceMethods
};

PyMODINIT_FUNC PyInit_surfaceModules(void){
    import_array();
    return PyModule_Create(&surfaceModules);
}

我用Python加载并使用这个模块，因此

import surfaceModules as sm

y=sm.squeezeSurfaces(surfaceVectors, magnitudes)

输出包括以下内容

f(89,75,84) = 64, v=(115,118,-63)
f(4,76,84) = 98, v=(-40,-123,-110)
f(5,76,84) = 101, v=(-63,85,61)
f(7,76,84) = 64, v=(115,118,-63)
f(12,76,84) = 98, v=(-40,-123,-110)
f(13,76,84) = 101, v=(-63,85,61)
f(15,76,84) = 64, v=(115,118,-63)
f(18,76,84) = 96, v=(0,0,0)
f(23,76,84) = 64, v=(-63,111,64)
f(26,76,84) = 96, v=(86,-63,85)
f(30,76,84) = 113, v=(0,0,0)
f(31,76,84) = 64, v=(-63,111,64)

但是，我在Python中获得了以下内容

surfaceVectors[84][76][31]
Out[20]: array([0., 0., 0.])

magnitudes[84][76][31]
Out[21]: 0.0

因此，C代码中PyArrayObject数组中的值与Python中相应numpy数组中的相应值不同。此外，范数值不是从中导出的3D向量所期望的值