Dec-05-2019, 11:55 PM
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.0So 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.