/*
 * create_circle.c - creates a circle of fixed radius
 *  The circle begins in the center of the image then
 *
 * MATLAB usage: y=create_shifted_circle(rows,cols,rad,rowcent,colcent)
 *
 */

#include <math.h>
#include "mex.h"
#define ROUND(A) ((fabs(A)) >= ((int)A+0.5) ? ((int)A+1.0) : ((int)A))  /* this is the rounding operation for + ints */

/* This is the original subroutine that performs the calculation...*/
double **create_shifted_circle(int rows, int cols, int rad, int rowcent, int colcent)
{
  int i,j;
  double D,**cir;
/*
  printf("in the C-function, rows = %d, cols = %d, rad = %d, rowcent=%d, colcent=%d\n",rows,cols,rad,rowcent,colcent);
*/  
  /* initialize output matrix */
  cir = (double **) mxCalloc(rows, sizeof(double *));
  cir[0]=(double *) mxCalloc(rows*cols, sizeof(double));
  for (i=1; i<rows; i++)
    cir[i]=cir[0] + i*cols;


  for (i=0;i<rows;i++)
    for (j=0;j<cols;j++)
       {
         D=sqrt((i-rowcent)*(i-rowcent)+(j-colcent)*(j-colcent));
		 
         if ((int)ROUND(D) == rad)
		    {
			 /*
			 printf("here, D=%f\n",D);
		     */
             cir[i][j]=1.0;
			}
       }
 return cir;
}

void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[])
{
  int rows,cols,rad,rowcent,colcent;
  int i,j,index=0;
  double **cir,**cir1,*y; /* array cir1 was created because MATLAB stores data in row ordered
						     form, so my C-result "cir" must be transposed before written out */

 /* Do some error checking */
 if (nrhs != 5)
   mexErrMsgTxt("Not enough input arguments.");
 else if (nlhs > 1)
   mexErrMsgTxt("Too many output arguments.");

/* Create a matrix for the return argument */
rows=mxGetScalar(prhs[0]);
cols=mxGetScalar(prhs[1]);
rad=mxGetScalar(prhs[2]);
rowcent=mxGetScalar(prhs[3])-1; /* C indices start at 0, vice 1 */
colcent=mxGetScalar(prhs[4])-1;
/*
printf("in the mex, rows = %d, cols = %d, rad = %d\n",rows,cols,rad);
*/

/* Assign pointers to the parameters */

cir=create_shifted_circle(rows,cols,rad,rowcent,colcent);

cir1 = (double **) mxCalloc(cols, sizeof(double *));
cir1[0]=(double *) mxCalloc(rows*cols, sizeof(double));
for (i=1; i<cols; i++)
    cir1[i]=cir1[0] + i*rows;
for (i=0;i<cols;i++)
  for (j=0;j<rows;j++)
	  cir1[i][j]=cir[j][i];

plhs[0]=mxCreateDoubleMatrix(rows,cols,mxREAL);  /* create output array */

y=(double *)mxGetPr(plhs[0]);

memcpy(y,cir1[0],rows*cols*sizeof(double));
}