Python Forum

I have a 4D numpy array, surfaceVectors, where the 4th dimension is the unnormalized vector (size and direction) at each voxel of a 3D volume. I make a corresponding 3D volume, magnitudes, where each element is the norm of the vector at the corresponding element in surfaceVectors. I make magnitudes as follows.

    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)

I then load a module, surfaceModules, which is a CInterface module, using

import surfaceModules as sm

y=sm.squeezeSurfaces(surfaceVectors, magnitudes)

squeezeSurfaces is a C function with the following form

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);
}

The output I get from

sm.squeezeSurfaces(surfaceVectors, magnitudes)

includes the following.

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)

However, I get the following in Python

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

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

So the values in the PyArrayObject arrays, in the C code, are not the same as the corresponding values in the corresponding numpy arrays in Python. Additionally, the norm values are not what would be expected given the 3D vectors from which they are derived.

I thought I should post the whole C code.

#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;
}

/* ==== Free a double *vector (vec of pointers) ========================== */
//void free_Carrayptrs(double **v) {
//free((char*) 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;
//                for (h=0; h<3; ++h)
//                    cout[i][j][k][h]= 1.0/h;
            }
	}  }

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

	return PyArray_Return(matout);

//	Py_INCREF(matin);
//	return PyArray_Return(matin);
}

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);
}