/***************************************************************************
  COPYRIGHT 1996 BIOMEDICAL COMPUTER LABORATORY, WASHINGTON UNIVERSITY, MO
****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <sys/types.h>
#ifndef WIN32
#include <sys/param.h>
#endif
#include "washu.h"
#include "fft3d.h"

#define DatOf1D(dat,x)                  *(dat + (x))
#define DatOf2D(dat,x,y,wid)            *(dat + (x) + (y)*wid)
#define MAXPATHLEN 256

extern char lognm[];
extern void WritePlane(int plane, float *outvol, float *inslice, int w, int h, int inx ,int iny);
extern int spline_(int *points, double *coord, float *samples, double *b_coef,
double *c_coef, double *d_coef);

/*************************************************************************************
  ROTNONE : Normal (optical sectioning) PSF volume creation.
	    This routine reads in the radial section of the PSF that was
	    generated by xcosm_psf and radially sweeps it 360 degrees
	    (centered at 0,0) by interpolating the values.

	    The program also checks for a symmetric psf, in which case
	    the above procedure need only be performed for N/2+1 rows	     
	    and then duplicated.

	    The resulting psf is centered at 0,0,0.	

 Author: Keith Doolittle
	 PSF algorithms developed by Jose-Angel Conchello
*************************************************************************************/
void rotnone(float *outimg, osm_ds *head)
{
char	time_str[MAXPATHLEN];
char	temp_str[MAXPATHLEN];
FILE	*fplog;
int   	oversmp,symmetric;
int   	N,NO2,Nxy,Nz,Nr,Nzwanted;
float 	deltaxy,deltar;
float 	*psfNN;
float  *outpsf, *rowdat;
int   	ix,iy,iz,ir,eta;
float 	r,x,y,fval;
float 	interp_val;
double  cord, sq_cord;
double  *bcf, *ccf, *dcf, *dr_cord;
time_t 	clock;

outpsf    = (float*)head->data;		/* XY slice data */
Nz 	  = head->ny;
Nr	  = head->nx;

oversmp   = head->xstart;		/* oversmp */
Nxy       = head->ystart;
symmetric = head->zstart;		/* symmetr */
deltar    = head->xlength;

N      = Nxy;			
NO2    = N/2;

/* Note: deltaxy and deltar are based on Nxy/2 = Nover2	*/

deltaxy = deltar*(float)oversmp;

if(symmetric == 0)
   Nzwanted = Nz;
else
   Nzwanted = Nz/2+1;

/* psf slice N x N 2-D plane */

if((psfNN=(float*)calloc(sizeof(float),N*N))==(float*)NULL) {
 fprintf(stderr,"rotnone; ERROR: Out of memory (6).\n");
 exit(1);
 }

rowdat = (float*) outpsf;

if((dr_cord=(double*)calloc(sizeof(double),Nr))==(double*)NULL) {
 fprintf(stderr,"rotnone; ERROR: Out of memory (6).\n");
 exit(1);
 }
if((bcf=(double*)calloc(sizeof(double),Nr))==(double*)NULL) {
 fprintf(stderr,"rotnone; ERROR: Out of memory (6).\n");
 exit(1);
 }
if((ccf=(double*)calloc(sizeof(double),Nr))==(double*)NULL) {
 fprintf(stderr,"rotnone; ERROR: Out of memory (6).\n");
 exit(1);
 }
if((dcf=(double*)calloc(sizeof(double),Nr))==(double*)NULL) {
 fprintf(stderr,"rotnone; ERROR: Out of memory (6).\n");
 exit(1);
 }

/* initialize radial coordinates array for the spline call */
 for(ir=0;ir<Nr;ir++) 
	dr_cord[ir] = ir * deltar;

for(iz=0;iz<Nzwanted;iz++) {

/* compute the cubic spline interpolation coefficients: bcf, ccf, dcf */
spline_(&Nr,dr_cord, rowdat,bcf,ccf,dcf);

 /* Sweep psf radial data into 1 plane */
 for(iy=0;iy<=NO2;iy++) {
  y = (float)iy * deltaxy;
  for(ix=0;ix<=NO2;ix++) {
   x = (float)ix * deltaxy;
   r = (float)sqrt(x*x + y*y);

   ir = (int) (r/deltar);
   cord = r - dr_cord[ir];
   sq_cord = cord * cord;


/* interpolated value based on cubic spline interpolation equation */
   interp_val = rowdat[ir] + bcf[ir]*cord + ccf[ir]*sq_cord + dcf[ir]*cord*sq_cord;

   /* Upper left */
   if((iy < NO2)&&(ix < NO2))
     DatOf2D(psfNN,ix,iy,N) = interp_val;
   /* Upper right */
   if(ix > 0)
     DatOf2D(psfNN,N-ix,iy,N) = interp_val;
   /* Lower left */
   if(iy > 0)
     DatOf2D(psfNN,ix,N-iy,N) = interp_val;
   /* Lower right */
   if((ix > 0)&&(iy > 0))
     DatOf2D(psfNN,N-ix,N-iy,N) = interp_val;
   }
  }

 /* Copy to 3-D PSF */

 WritePlane(iz,outimg,psfNN,N,N,N,N);
 if((symmetric != 0)&&(iz > 0)&&(iz < Nzwanted-1))
	WritePlane(Nz-iz,outimg,psfNN,N,N,N,N);

 /* Increment row counter to point to next row */
 rowdat += Nr;

 if((fplog=fopen(lognm,"a"))==(FILE*)NULL)
   fprintf(stderr,"rotnone; WARNING: Can't open logfile `%s' for write/append (rotnormal)\n",lognm);
 else {
      clock = time(NULL);
      strcpy(time_str,asctime(localtime(&clock)));
      time_str[strlen(time_str)-1] = '\0';
      sprintf(temp_str," Completed plane %d",iz);
      strcat(time_str,temp_str);
      fprintf(fplog,"rotnone; %s\n",time_str);
      fclose(fplog);
      }
 }

}

