/*
 * $Id: spm_resels_vol.c 4453 2011-09-02 10:47:25Z guillaume $
 * John Ashburner
 */
 
/*
See Worsley et al. (1996), Human Brain Mapping 4:58-73 for a description
of what it does.
*/

#include <math.h>
#include "mex.h"
#include "spm_mapping.h"

static void resel_fun(int *curr, int *prev, /* current and previous planes */
    int m, int n, /* image dimensions */
    int *P,   /* # points */
    int E[3], /* # edges  */
    int F[4], /* # faces  */
    int *C)   /* # cubes  */
{
    int p=0,ex=0,ey=0,ez=0,fxy=0,fxz=0,fyz=0,c=0;
    int i, j;
    for(i=1; i<m; i++)
        for(j=1; j<n; j++)
        {
            int o = i+j*m;

            if (curr[o])
            {
                p ++;

                /* A simple way of doing it... 
                if (curr[o-1]) ex ++;
                if (curr[o-m]) ey ++;
                if (prev[o  ]) ez ++;
                if (curr[o-1] && curr[o-m] && curr[o-m-1]) fxy ++;
                if (curr[o-1] && prev[o  ] && prev[o  -1]) fxz ++;
                if (curr[o-m] && prev[o  ] && prev[o  -m]) fyz ++;
                if (curr[o-1] && curr[o-m] && curr[o-m-1] && prev[o]
                 && prev[o-1] && prev[o-m] && prev[o-m-1]) c ++;
                */

                /* It could be optimized much further with "if else" statements
                   but it was beginning to hurt my head */
                if (curr[o-1])
                {
                    ex ++;
                    if (curr[o-m] && curr[o-m-1])
                    {
                        fxy ++;
                        if (prev[o] && prev[o-1] && prev[o-m] && prev[o-m-1]) c ++;
                    }
                }
                if (curr[o-m])
                {
                    ey ++;
                    if (prev[o  ] && prev[o  -m]) fyz ++;
                }
                if (prev[o  ])
                {
                    ez ++;
                    if (curr[o-1] && prev[o  -1]) fxz ++;
                }
            }
        }
    *P   += p;
    E[0] += ex;
    E[1] += ey;
    E[2] += ez;
    F[0] += fxy;
    F[1] += fxz;
    F[2] += fyz;
    *C   += c;
}


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    MAPTYPE *map, *get_maps();
    mwSize m,n,k;
    mwIndex i;
    int nn, E[3], F[3], P, C;
    double *R, r[3], *img;
    int *curr, *prev, *tmpp;

    if (nrhs != 2 || nlhs > 1)
    {
        mexErrMsgTxt("Incorrect usage.");
    }

    map = get_maps(prhs[0], &nn);
    if (nn!=1)
    {
        free_maps(map, nn);
        mexErrMsgTxt("Bad image handle dimensions.");
    }

    if (!mxIsNumeric(prhs[1]) || mxIsComplex(prhs[1]) || !mxIsDouble(prhs[1]))
    {
        free_maps(map, 1);
        mexErrMsgTxt("Second argument must be numeric, real, full and double.");
    }

    if (mxGetM(prhs[1])*mxGetN(prhs[1]) != 3)
    {
        free_maps(map, nn);
        mexErrMsgTxt("Second argument must contain three elements.");
    }

    r[0] = 1.0/mxGetPr(prhs[1])[0];
    r[1] = 1.0/mxGetPr(prhs[1])[1];
    r[2] = 1.0/mxGetPr(prhs[1])[2];

    plhs[0] = mxCreateDoubleMatrix(4,1,mxREAL);
    R       = mxGetPr(plhs[0]);

    m = map->dim[0];
    n = map->dim[1];
    k = map->dim[2];

    curr = (int *)mxCalloc((m+1)*(n+1),sizeof(int)); /* current plane */
    prev = (int *)mxCalloc((m+1)*(n+1),sizeof(int)); /* previous plane */
    img  = (double *)mxCalloc((m+1)*(n+1),sizeof(double));

    P = C = E[0] = E[1] = E[2] = F[0] = F[1] = F[2] = 0;
    for(i=0; i<k; i++)
    {
        int i1, j1;
        static double mat[16] = {
            1.0, 0.0, 0.0, 0.0,
            0.0, 1.0, 0.0, 0.0,
            0.0, 0.0, 1.0, 0.0,
            0.0, 0.0, 0.0, 1.0};
        mat[14] = i+1.0;
        slice(mat, img, m, n, map, 0, 0);
        for(i1=0;i1<m; i1++)
            for(j1=0; j1<n; j1++)
                if (mxIsFinite(img[i1+m*j1]) && img[i1+m*j1])
                    curr[i1+1+(m+1)*(j1+1)] = 1;
                else
                    curr[i1+1+(m+1)*(j1+1)] = 0;

        /* count edges, faces etc */
        resel_fun(curr,prev,m+1,n+1, &P, E, F, &C);

        /* make current plane previous */
        tmpp = prev; prev = curr; curr = tmpp;
    }

    (void)mxFree((char *)curr);
    (void)mxFree((char *)prev);
    (void)mxFree((char *)img);
    free_maps(map, 1);

    R[0] = P - (E[0]+E[1]+E[2])+(F[0]+F[1]+F[2])-C;
    R[1] = (E[0]-F[0]-F[1]+C)*r[0] + (E[1]-F[0]-F[2]+C)*r[1] + (E[2]-F[1]-F[2]+C)*r[2];
    R[2] = (F[0]-C)*r[0]*r[1] + (F[1]-C)*r[0]*r[2] + (F[2]-C)*r[1]*r[2];
    R[3] = C*r[0]*r[1]*r[2];
}