float *x1, *x2, **ys;
float *yy1, *ys1, *YS1, *cys1;
void magnify1(float *image, int rows, int cols)
{
float result;
float frow, fcol;
int i, r, c;
//calloc for 1 D
for (row = 0; row < rows; row++)
x1[row] = (float) row;
for (col = 0; col < cols; col++)
{
for (row = 0; row < rows; row++)
yy1[row] = image[(row + 37)*lenslets_size + col + 37]; // for each row you have an array of the col number and an array of pixel value
spline( x1,yy1,rows,0,0, ys1); //make the ys1 array
memcpy(YS1+col*reduced_size, ys1, sizeof(float)*reduced_size);
// printf("here \n");
// for ( i = 0; i < 10; i ++)
// printf("here %f \n", *(YS1+col*reduced_size + i) );
}
//return (YS1);
}
void getmemory(int rows, int cols)
{
if ((x1 = (float *)calloc(reduced_size,sizeof(float))) == NULL)
{
printf("no x1 memory\n");
exit(0);
}
if ((yy1 = (float *)calloc(reduced_size,sizeof(float))) == NULL)
{
printf("no x2 memory\n");
exit(0);
}
if ((ys1 = (float *)calloc(reduced_size ,sizeof(float))) == NULL)
{
printf("no ys1 memory\n");
exit(0);
}
if ((YS1 = (float *)calloc(reduced_size *reduced_size ,sizeof(float))) == NULL)
{
printf("no YS1 memory\n");
exit(0);
}
if ((cys1 = (float *)calloc(waves*reduced_size*reduced_size ,sizeof(float))) == NULL)
{
printf("no cys1 memory\n");
exit(0);
}
}
void spline(float x[], float y[], int n, float yp1, float ypn, float y2[])
{
int i,k;
float p,qn,sig,un,*u;
if ( !( u = (float *)calloc (2*n,sizeof (float)) ))
printf( "calloc failed\n");
if (yp1 > 0.99e30) //The lower boundary condition is set either to be "natural"
y2[0]=u[0]=0.0;
else { //or else to have a specified first derivative.
y2[0] = -0.5;
u[0]=(3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
}
for (i=1;i<=n-2;i++) { //This is the decomposition loop of the tridiagonal algorithm. y2 and u are used for temporary
//storage of the decomposed factors.
sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
p=sig*y2[i-1]+2.0;
y2[i]=(sig-1.0)/p;
u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
}
n = n-1;
if (ypn > 0.99e30) //The upper boundary condition is set either to be "natural"
qn=un=0.0; //or else to have a specified first derivative.
else {
qn=0.5;
un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
}
y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
for (k=n;k>=0;k--) //This is the backsubstitution loop of the tridiagonal algorithm.
y2[k]=y2[k]*y2[k+1]+u[k];
// printf("here\n");
free (u);
}