/*This is a mex file for caculate the regularization term. \int_{B(x;wsize)(1-u_k(x))};
input U is a 4-D array which stands for u(x,y,z,k). 
wsize/2 is the radius of a ball.
The output is a 4-D array represents the regularization function R(x,y,z,k).
This code needs the CImg Libaray to support.
Author: Jun Liu. Dec. 2010.
*/ 
#include <mex.h>
#define cimg_display 0
#define cimg_plugin "cimgmatlab.h"
#include "CImg.h"
//#define cimg_plugin "/media/Study/C&C++/CImg/CImg-1.3.4/CImg-1.3.4/plugins/cimgmatlab.h"
//#include "/media/Study/C&C++/CImg/CImg-1.3.4/CImg-1.3.4/CImg.h"
//#define cimg_plugin "/media/Study/C&C++/CImg/CImg-1.3.4/CImg-1.3.4/plugins/cimgmatlab.h"
//#include"/media/Study/C&C++/CImg/CImg-1.3.4/CImg-1.3.4/CImg.h"
using namespace cimg_library;
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{

       CImg<> U(prhs[0],false);
       CImg<> pfVecParameters(prhs[1],true);
       int wsize= (int) pfVecParameters[0];
       const int wsizeleft= -wsize/2,wsizeright=(wsize+1)/2;
       const int Nx=U.width(),Ny=U.height(),Nz=U.depth(),Nc=U.spectrum();
       CImg<> R(Nx,Ny,Nz,Nc);
       cimg_forXYZC(U,x,y,z,c)
       	{
	       float dist=0.0f;
		   for (int iz=wsizeleft;iz<wsizeright;iz++ )
		       for (int iy=wsizeleft;iy<wsizeright;iy++ )
			  for (int ix=wsizeleft;ix<wsizeright;ix++ )
                              if ((float)(ix*ix+iy*iy+iz*iz)<(float)(wsize*wsize)/4.0f)
				   {
                                         int Ix=x+ix,Iy=y+iy,Iz=z+iz;
					 if (Ix<0) Ix=-Ix; else if (Ix>Nx-1) Ix= 2*(Nx-1)-Ix;
					 if (Iy<0) Iy=-Iy; else if (Iy>Ny-1) Iy= 2*(Ny-1)-Iy;
				         if (Iz<0) Iz=-Iz; else if (Iz>Nz-1) Iz= 2*(Nz-1)-Iz;
			                  dist+=(1.0f-U(Ix,Iy,Iz,c));

				   }
                  R(x,y,z,c)=dist;
           }

   plhs[0]=R.toMatlab();
}
