Page MenuHomec4science
Files F122453851

spm_get_lm.c
No OneTemporary
Actions

File Metadata

Created
Thu, Jul 17, 23:16

spm_get_lm.c
View Options

/*
* $Id: spm_get_lm.c 4453 2011-09-02 10:47:25Z guillaume $
* Jesper Andersson
*/
/****************************************************************
**
** Routine that identifies which voxels in a list of coordinates
** that are local maxima, and returns a list of indicies into
** the coordinate list for those maxima.
**
***************************************************************/
#include <math.h>
#include <limits.h>
#include <string.h>
#include "mex.h"
/* Macros */
#ifndef MIN
#define MIN(A,B) ((A) > (B) ? (B) : (A))
#endif
#ifndef MAX
#define MAX(A,B) ((A) > (B) ? (A) : (B))
#endif
/* get_index */
mwSignedIndex get_index(mwIndex x, mwIndex y, mwIndex z, mwSize dim[3])
{
if (x < 1 || x > dim[0] || y < 1 || y > dim[1] || z < 1 || z > dim[2])
{
return(-1);
}
else
{
return((z-1)*dim[0]*dim[1]+(y-1)*dim[0]+x-1);
}
}
/* is_maxima */
int is_maxima(double *v,
mwSize dim[3],
mwIndex x,
mwIndex y,
mwIndex z,
unsigned int cc)
{
mwSignedIndex ii = 0, i = 0;
double cv = 0.0;
if ((ii=get_index(x,y,z,dim))<0) {return(0);}
cv = v[ii];
if (cc >= 6)
{
if ((i=get_index(x+1,y,z,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x-1,y,z,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x,y+1,z,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x,y-1,z,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x,y,z+1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x,y,z-1,dim))>0 && v[i] > cv) {return(0);}
}
if (cc >= 18)
{
if ((i=get_index(x+1,y+1,z,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x+1,y-1,z,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x+1,y,z+1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x+1,y,z-1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x-1,y+1,z,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x-1,y-1,z,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x-1,y,z+1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x-1,y,z-1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x,y+1,z+1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x,y+1,z-1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x,y-1,z+1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x,y-1,z-1,dim))>0 && v[i] > cv) {return(0);}
}
if (cc == 26)
{
if ((i=get_index(x+1,y+1,z+1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x+1,y+1,z-1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x+1,y-1,z+1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x+1,y-1,z-1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x-1,y+1,z+1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x-1,y+1,z-1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x-1,y-1,z+1,dim))>0 && v[i] > cv) {return(0);}
if ((i=get_index(x-1,y-1,z-1,dim))>0 && v[i] > cv) {return(0);}
}
return(1);
}
/* get_maxima */
unsigned int get_maxima(double *vol,
mwSize vdim[3],
double *list,
mwSize nlist,
unsigned int cc,
mwIndex **ldx)
{
mwIndex i = 0, j = 0;
mwIndex ix = 0, iy = 0, iz = 0;
mwSize ldx_sz = 0;
mwIndex ldx_n = 0;
*ldx = (mwIndex *) mxCalloc((ldx_sz = 1024),sizeof(mwIndex));
for (i=0, j=0; i<nlist; i++, j+=3)
{
/*
** Casting of double to int isn't properly defined in C
** (i.e. wether it results in truncation or rounding),
** hence I add a small offset (0.1) to make sure it
** works either way.
*/
ix = ((mwIndex) (list[j]+0.1));
iy = ((mwIndex) (list[j+1]+0.1));
iz = ((mwIndex) (list[j+2]+0.1));
if (get_index(ix,iy,iz,vdim) > 0)
{
if (is_maxima(vol,vdim,ix,iy,iz,cc))
{
if (ldx_n >= ldx_sz)
{
*ldx = (mwIndex *) mxRealloc(*ldx,(ldx_sz += 1024)*sizeof(mwIndex));
}
(*ldx)[ldx_n] = i+1;
ldx_n++;
}
}
}
return(ldx_n);
}
/* Gateway function */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mwIndex i = 0, j = 0, k = 0;
int tmpint = 0;
const mwSize *pdim = NULL;
mwSize ndim = 0;
mwSize vdim[3];
mwSize ln = 0, lm = 0;
mwSize n_lindex = 0;
mwIndex *lindex = NULL;
double *vol = NULL;
double *lp = NULL;
double *list = NULL;
double *plindex = NULL;
if (nrhs < 2) mexErrMsgTxt("Not enough input arguments.");
if (nrhs > 2) mexErrMsgTxt("Too many input arguments.");
if (nlhs < 1) mexErrMsgTxt("Not enough output arguments");
if (nlhs > 1) mexErrMsgTxt("Too many output arguments.");
/* Get binary map. */
if (!mxIsNumeric(prhs[0]) || mxIsComplex(prhs[0]) || mxIsSparse(prhs[0]) || !mxIsDouble(prhs[0]))
{
mexErrMsgTxt("spm_get_lm: VOL must be numeric, real, full and double");
}
ndim = mxGetNumberOfDimensions(prhs[0]);
if (ndim != 3 && ndim != 2)
{
mexErrMsgTxt("spm_get_lm: VOL must 2- or 3-dimensional");
}
pdim = mxGetDimensions(prhs[0]);
vdim[0] = pdim[0]; vdim[1] = pdim[1];
if (ndim == 2) { vdim[2] = 1; } else { vdim[2] = pdim[2]; }
vol = mxGetPr(prhs[0]);
/* Get list of coordinates */
if (!mxIsNumeric(prhs[1]) || mxIsComplex(prhs[1]) || mxIsSparse(prhs[1]) || !mxIsDouble(prhs[1]))
{
mexErrMsgTxt("spm_get_lm: L must be numeric, real, full and double");
}
lm = mxGetM(prhs[1]);
ln = mxGetN(prhs[1]);
if (!((lm==3 && ndim==3) || (lm==3 && ndim==2) || (lm==2 && ndim==2)))
{
mexErrMsgTxt("spm_get_lm: L must be 3xn (or 2xn) list of voxel coordinates");
}
lp = mxGetPr(prhs[1]);
if (lm==3 && ndim==2) /* Make sure all z-coordinates equal 1 if 3xn list used with 2D map. */
{
for (i=0, j=0; i<ln; i++, j+=3)
{
tmpint = ((int) lp[j+2]+0.1);
if (tmpint != 1)
{
mexErrMsgTxt("spm_get_lm: z-coordinate must be 1 when using 3xn list with 2D map");
}
}
}
list = (double *) mxCalloc(3*ln,sizeof(double));
if (lm==2) /* Extend to 3xn list if 2xn list was supplied with 2D-map. */
{
for (i=0, j=0, k=0; i<ln; i++, j+=3, k+=2)
{
list[j] = lp[k]; list[j+1] = lp[k+1]; list[j+2] = 1.0;
}
}
else
{
memcpy(list,lp,3*ln*sizeof(double));
}
/* Find list if indicies to local maxima. */
n_lindex = get_maxima(vol,vdim,list,ln,18,&lindex);
/* Turn indicies into doubles in a MATLAB array. */
plhs[0] = mxCreateDoubleMatrix(1,n_lindex,mxREAL);
plindex = mxGetPr(plhs[0]);
for (i=0; i<n_lindex; i++) {plindex[i] = ((double) lindex[i]);}
mxFree(list);
mxFree(lindex);
}

Event Timeline

c4science · Help