text: /*
   * $Id: spm_project.c 4453 2011-09-02 10:47:25Z guillaume $
   * John Ashburner
   */
  
  #include <math.h>
  #include <stdio.h>
  #include "mex.h"
  #define RINT(A) floor((A)+0.5)
  #define MAX(A, B)   ((A) > (B) ? (A) : (B))
  #define MIN(A, B)   ((A) < (B) ? (A) : (B))
  
  void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
  {
      double *spm,*l,*v,*dim;
      mwSize m,n;
      int    m1,i,j,k, o;
      int    x,y,z,xdim,ydim,zdim;
      double q;
      double *DXYZ, *CXYZ;
      int    DX, DY, DZ, CX, CY, CZ;
  
      if ( !(nrhs == 3 || nrhs == 5) || nlhs > 1) mexErrMsgTxt("Incorrect usage.");
  
      for(k=0; k<nrhs; k++)
          if (!mxIsNumeric(prhs[k]) || mxIsComplex(prhs[k]) ||
              mxIsSparse(prhs[k]) || !mxIsDouble(prhs[k]))
              mexErrMsgTxt("Arguments must be numeric, real, full and double.");
  
      /* The values */
      n    = mxGetNumberOfElements(prhs[0]);
      v    = mxGetPr(prhs[0]);
  
      /* The co-ordinates */
      if ((mxGetN(prhs[1]) != n) || (mxGetM(prhs[1]) != 3))
          mexErrMsgTxt("Incompatible size for locations matrix.");
      l    = mxGetPr(prhs[1]);
  
      /* Dimension information */
      if (mxGetNumberOfElements(prhs[2]) != 5)
          mexErrMsgTxt("Incompatible size for dimensions vector.");
      dim  = mxGetPr(prhs[2]);
      xdim = (int) (fabs(dim[0]) + 0.99);
      ydim = (int) (fabs(dim[1]) + 0.99);
      zdim = (int) (fabs(dim[2]) + 0.99);
      m    = (int) (dim[3]);
      m1   = (int) (dim[4]);
  
      if (nrhs == 3)
      {
        DX = 182;
        DY = 218;
        DZ = 182;
        CX =  89;
        CY = 125;
        CZ =  71;
      }
      else
      {
        if (mxGetNumberOfElements(prhs[3]) != 3)
          mexErrMsgTxt("Incompatible size for DXYZ vector.");
        DXYZ = mxGetPr(prhs[3]);
        DX = (int) (DXYZ[0]);
        DY = (int) (DXYZ[1]);
        DZ = (int) (DXYZ[2]);
        if (mxGetNumberOfElements(prhs[4]) != 3)
          mexErrMsgTxt("Incompatible size for CXYZ vector.");
        CXYZ = mxGetPr(prhs[4]);
        CX = (int) (CXYZ[0]);
        CY = (int) (CXYZ[1]);
        CZ = (int) (CXYZ[2]);
      }
  
      plhs[0] = mxCreateDoubleMatrix(m,(mwSize)m1,mxREAL);
      spm     = mxGetPr(plhs[0]);
  
      if (m == DY+DX && m1 == DZ+DX) /* MNI Space */
      {
          /* go though point list */
          for (i = 0; i < n; i++)
          {
              x = (int)RINT(l[i*3 + 0]) + CX;
              y = (int)RINT(l[i*3 + 1]) + CY;
              z = (int)RINT(l[i*3 + 2]) + CZ;
  
              if (2*CX-x-xdim/2>=0 && 2*CX-x+xdim/2<DX && y-ydim/2>=0 && y+ydim/2<DY) /* transverse */
              {
                  q = v[i];
                  for (j = -ydim/2; j <= ydim/2; j++)
                      for (k = -xdim/2; k <= xdim/2; k++)
                      {
                          o = j + y + (k + 2*CX-x)*m;
                          if (spm[o]<q) spm[o] = q;
                      }
              }
  
              if (z-zdim/2>=0 && z+zdim/2<DZ && y-ydim/2>=0 && y+ydim/2<DY) /* sagittal */
              {
                  q = v[i];
                  for (j = -ydim/2; j <= ydim/2; j++)
                      for (k = -zdim/2; k <= zdim/2; k++)
                      {
                          o = j + y + (DX + k + z)*m;
                          if (spm[o]<q) spm[o] = q;
                      }
              }
  
              if (x-xdim/2>=0 && x+xdim/2<DX && z-zdim/2>=0 && z+zdim/2<DZ) /* coronal */
              {
                  q = v[i];
                  for (j = -xdim/2; j <= xdim/2; j++)
                      for (k = -zdim/2; k <= zdim/2; k++)
                      {
                          o = DY + j + x + (DX + k + z)*m;
                          if (spm[o]<q) spm[o] = q;
                      }
              }
          }
      }
      else if (m == 360 && m1 == 352) /* old code for the old MIP matrix */
      {
      for (i = 0; i < n; i++) {
          x = (int) l[i*3 + 0];
          y = (int) l[i*3 + 1];
          z = (int) l[i*3 + 2];
      
          /* transverse */
          q = MAX(v[i], spm[(124 + y) + (104 - x)*m]);
          for (j = 0; j < ydim; j++) {
              for (k = 0; k < xdim; k++) {
                  spm[124 + j + y + (104 + k - x)*m] = q;
              }
          }
      
          /* sagittal */
          q = MAX(v[i], spm[(124 + y) + (240 + z)*m]);
          for (j = 0; j < ydim; j++) {
              for (k = 0; k < zdim; k++) {
                  spm[124 + j + y + (238 + k + z)*m] = q;
              }
          }
      
          /* coronal */
          q = MAX(v[i], spm[(276 + x) + (240 + z)*m]);
          for (j = 0; j < xdim; j++) {
              for (k = 0; k < zdim; k++) {
                  spm[276 + j + x + (238 + k + z)*m] = q;
              }
          }
      }
      }
      else 
          mexErrMsgTxt("Wrong sized MIP matrix");
  }
  