#ifndef lint static char sccsid[]="@(#)spm_clusters.c 2.2 JB Poline 99/03/03"; #endif /* spm_plateau.c JB Poline 10/11/94 % Return cluster index for a point list % FORMAT [A] = spm_cluster2(L) % L - locations [x y x]' {in vox} % % A - cluster index or region number %____________________________________________________________________________ % % spm_clusters characterizes a point list of voxel values (X) and their % locations (L) in terms of edge, face and vertex connected subsets, returning a % list of indices in A, such that X(i) belongs to cluster A(i) {using an 18 % connectivity scheme) % */ #include #include "mex.h" /* Input Output Arguments */ #define XYZ prhs[0] #define IND plhs[0] #include "connex.h" /* where typedef struct cc_position_3d_str { float t, sum, label, *ptr_flt; int x, y, z, dx, dy, dz, size; } cc_position_3d; static cc_position_3d pos3; is defined */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double *L; float *vol, fl_tmp; int i,j,d,r; int max_x, max_y, max_z, min_x, min_y, min_z, dx, dy, dz, dxdy, *ind, index; int *x,*y,*z; if ((nrhs != 1) || (nlhs > 1)) mexErrMsgTxt("Inappropriate usage."); /* Assign pointers to the parameters */ L = mxGetPr(XYZ); i = mxGetM(XYZ); if(i != 3) mexErrMsgTxt(" XYZ M-dimension should be 3 "); r = mxGetN(XYZ); if(!r) mexErrMsgTxt("no points in the volume in spm_max"); /* set up point list in x y z */ x = (int *) mxCalloc (r,sizeof( int )); y = (int *) mxCalloc (r,sizeof( int )); z = (int *) mxCalloc (r,sizeof( int )); ind = (int *) mxCalloc (r,sizeof( int )); if(!x || !y || !z || !ind) mexErrMsgTxt("\n memory alloc pb in loc_max "); max_x = (int) floor(L[0] + 0.5); max_y = (int) floor(L[1] + 0.5); max_z = (int) floor(L[2] + 0.5); min_x = (int) floor(L[0] + 0.5); min_y = (int) floor(L[1] + 0.5); min_z = (int) floor(L[2] + 0.5); d = 0; for (i = 0; i < 3*r; i = i + 3) { x[d] = (int) floor(L[i + 0] + 0.5); if (x[d] > max_x) max_x = x[d]; if (x[d] < min_x) min_x = x[d]; y[d] = (int) floor(L[i + 1] + 0.5); if (y[d] > max_y) max_y = y[d]; if (y[d] < min_y) min_y = y[d]; z[d] = (int) floor(L[i + 2] + 0.5); if (z[d] > max_z) max_z = z[d]; if (z[d] < min_z) min_z = z[d]; d++; } /* create a volume of the points plus a bounding box */ dx = max_x - min_x + 3; dy = max_y - min_y + 3; dz = max_z - min_z + 3; dxdy = dx*dy; vol = (float *) mxCalloc((dx*dy*dz), sizeof(*vol)); if(!vol) mexErrMsgTxt("\n memory alloc pb in loc_max "); /* put the points in the volume */ for( i=0; i