#include <time.h>

#include <math.h>
#include "splat.h"

float alphaInt,alphaDeriv;
char normalEstimationMethod[100];
extern void Metrics(int, double *);

float int_filter[6];
float der_filter[6];
int   int_filter_extent;
int   der_filter_extent;

/***********************************************************************/
void NormalizeV(Vector *v) {
  double mag;

  mag=sqrt(Sqr(v->x)+Sqr(v->y)+Sqr(v->z));
  if(mag>TINY) {
    v->x/=mag; v->y/=mag; v->z/=mag;
  }
}

void InitMat3(Matrix3 mat) {
  int i,j;

  for(i=0;i<3;i++)
    for(j=0;j<3;j++)
      mat[i][j]=0.;
}

void SetRotationMatrix(Matrix3 mat,double angle,char axis) {
  InitMat3(mat);
  switch(axis) {
    case 'x' :
      mat[1][1]=(float)cos((double)angle); 
      mat[1][2]=-(float)sin((double)angle); 
      mat[2][1]=(float)sin((double)angle); 
      mat[2][2]=(float)cos((double)angle); 
      mat[0][0]=1.0;
    break;
    case 'y' :
      mat[0][0]=(float)cos((double)angle); 
      mat[0][2]=-(float)sin((double)angle); 
      mat[2][0]=(float)sin((double)angle); 
      mat[2][2]=(float)cos((double)angle); 
      mat[1][1]=1.0;
    break;
    case 'z' :
      mat[0][0]=(float)cos((double)angle); 
      mat[1][0]=(float)sin((double)angle); 
      mat[0][1]=-(float)sin((double)angle); 
      mat[1][1]=(float)cos((double)angle); 
      mat[2][2]=1.0;
    break;
  }
}

void Mat3Mult(Matrix3 mat1,Matrix3 mat2,Matrix3 matout) {
  int i,j,k;

  for(i=0;i<3;i++)
  {
    for(j=0;j<3;j++)
    {
      matout[i][j]=0;
      for(k=0;k<3;k++)
        matout[i][j]+=mat1[i][k]*mat2[k][j];
    }
  }
}

Vector VecMat3Mult(Vector *v,Matrix3 m) {
   Vector v1;
   
   v1.x=m[0][0]*v->x+m[1][0]*v->y+m[2][0]*v->z;
   v1.y=m[0][1]*v->x+m[1][1]*v->y+m[2][1]*v->z;
   v1.z=m[0][2]*v->x+m[1][2]*v->y+m[2][2]*v->z;

   return(v1);
}      

double Dot(Vector *v1,Vector *v2) {
  return(v1->x*v2->x+v1->y*v2->y+v1->z*v2->z);
}

Vector Cross(Vector *v1,Vector *v2) {
  Vector outVec;

  outVec.x=v1->y*v2->z-v1->z*v2->y;
  outVec.y=v1->z*v2->x-v1->x*v2->z;
  outVec.z=v1->x*v2->y-v1->y*v2->x;
  return(outVec);
}  

Vector ProjectVector(Vector *v,Vector *n) {
  double dot;
  Vector vp;

  dot=Dot(v,n);
  vp.x=v->x-n->x*dot;
  vp.y=v->y-n->y*dot;
  vp.z=v->z-n->z*dot;
  
  return(vp);
}    
 
void ComputeViewPlaneVectors(ViewParameters *vpars) {
  double eyeDist,denom;
  Matrix3 mat;
  Vector VPU,VPV,VPN,vec2;
  int i,j;

  vpars->VPV=ProjectVector(&vpars->VUP,&vpars->VPN);  
  NormalizeV(&vpars->VPV);

  vpars->VPU=Cross(&vpars->VPN,&vpars->VPV);
  NormalizeV(&vpars->VPU);

  if(strstr(vpars->splatMethod,"splat")) {
    VPU=vpars->VPU;
    VPV=vpars->VPV;
    VPN=vpars->VPN;

    denom=VPN.x*(VPU.z*VPV.y-VPU.y*VPV.z)+VPN.y*(VPU.x*VPV.z-VPU.z*VPV.x)+
      VPN.z*(VPU.y*VPV.x-VPU.x*VPV.y);

    mat[0][0]=(VPV.z*VPU.y-VPV.y*VPU.z)/denom;
    mat[0][1]=(VPN.y*VPV.z-VPN.z*VPV.y)/denom;
    mat[0][2]=(VPU.y*VPN.z-VPU.z*VPN.y)/denom;
    mat[1][0]=(VPU.z*VPV.x-VPU.x*VPV.z)/denom;
    mat[1][1]=(VPN.z*VPV.x-VPV.z*VPN.x)/denom;
    mat[1][2]=(VPN.x*VPU.z-VPU.x*VPN.z)/denom;
    mat[2][0]=(VPV.y*VPU.x-VPV.x*VPU.y)/denom;
    mat[2][1]=(VPV.y*VPN.x-VPN.y*VPV.x)/denom;
    mat[2][2]=(VPN.y*VPU.x-VPU.y*VPN.x)/denom;

    for(i=0;i<3;i++)
      for(j=0;j<3;j++)
        vpars->volRotationMat[i][j]=mat[i][j];  
  }

  if(vpars->viewAngle!=0) {
    eyeDist=(vpars->imgHeight/2.0)/tan(Rad((abs(vpars->viewAngle)/2.0)));
    vpars->EYE.x=vpars->VRP.x-eyeDist*vpars->VPN.x;
    vpars->EYE.y=vpars->VRP.y-eyeDist*vpars->VPN.y;
    vpars->EYE.z=vpars->VRP.z-eyeDist*vpars->VPN.z;
  }
}

int CreateGradientVolumes(Volume *vol,Volume *dxVol,Volume *dyVol,
			  Volume *dzVol,Volume *magVol,ViewParameters *vpars) {
  int x,y,z,n,nxny,nx,ny,nz,nx1,ny1,nz1,i;
  float *volPtr,*magData;
  float *dxData,*dyData,*dzData;
  double xGrad,yGrad,zGrad,mag;
  FILE *fpdx,*fpdy,*fpdz,*fpmag,*fp;

  nx=vol->nx;  ny=vol->ny;  nz=vol->nz;
  nx1=nx-1; ny1=ny-1; nz1=nz-1;
  nxny=nx*ny;  
  n=nx*ny*nz;

  dxData=(float*)calloc(nx*ny*nz,sizeof(float));
  dyData=(float*)calloc(nx*ny*nz,sizeof(float));
  dzData=(float*)calloc(nx*ny*nz,sizeof(float));  
  magData=(float*)calloc(nx*ny*nz,sizeof(float));  

  dxVol->fdata=dxData;
  dyVol->fdata=dyData;
  dzVol->fdata=dzData;
  magVol->fdata=magData;

  dxVol->nx=nx; dyVol->nx=nx; dzVol->nx=nx; magVol->nx=nx;
  dxVol->ny=ny; dyVol->ny=ny; dzVol->ny=ny; magVol->nx=ny;
  dxVol->nz=nz; dyVol->nz=nz; dzVol->nz=nz; magVol->nx=nz;

  if(!strcmp(vpars->gradientVolumeFlag,"use")) {
    if(!(fpdx =fopen("gradVolDx","r"))|| 
       !(fpdy =fopen("gradVolDy","r"))|| 
       !(fpdz =fopen("gradVolDz","r"))||
       !(fpmag=fopen("gradVolMag","r")))     {
      printf("\nCannot open gradient volumes\n");
      return(0);
    }
    if((fread(dxData,sizeof(float),n,fpdx)<n) ||
       (fread(dyData,sizeof(float),n,fpdy)<n) ||
       (fread(dzData,sizeof(float),n,fpdz)<n) ||
       (fread(magData,sizeof(float),n,fpmag)<n))   {
      printf("\nGradient volume size is different from volume size\n");
      return(0);
    }
    fclose(fpdx); fclose(fpdy); fclose(fpdz); fclose(fpmag);
  }    
  else {
    if(!strcmp(vpars->derivMethod,"central difference")) {
      volPtr=vol->fdata+nxny+ny+1;
      dxData+=nxny+ny+1;  dyData+=nxny+ny+1;  dzData+=nxny+ny+1;  
      magData+=nxny+ny+1;  
      for(z=1;z<nz1;z++)     {
        for(y=1;y<ny1;y++)   {
          for(x=1;x<nx1;x++) {  
            xGrad=(*(volPtr-1)-*(volPtr+1))/2.;        
            yGrad=(*(volPtr-nx)-*(volPtr+nx))/2.;      
            zGrad=(*(volPtr-nxny)-*(volPtr+nxny))/2.;  
            mag=sqrt(Sqr(xGrad)+Sqr(yGrad)+Sqr(zGrad));
            *magData=mag;
            *dxData=xGrad; 
            *dyData=yGrad; 
            *dzData=zGrad; 
             volPtr++; dxData++; dyData++; dzData++; magData++;
          }
          volPtr+=2; dxData+=2; dyData+=2; dzData+=2;
        }    
        volPtr+=2*nx; dxData+=2*nx; dyData+=2*nx; dzData+=2*nx;
      }
    }
    else   {
      printf("\nDon't know how to make %s filter\n",vpars->derivMethod);
      return(0);
    }
  }
        
  if(!strcmp(vpars->gradientVolumeFlag,"create"))  {
    fpdx=fopen("gradVolDx","w"); 
    fpdy=fopen("gradVolDy","w"); 
    fpdz=fopen("gradVolDz","w");
    fpmag=fopen("gradVolMag","w");

    fwrite(dxVol->fdata,sizeof(float),n,fpdx);
    fwrite(dyVol->fdata,sizeof(float),n,fpdy);
    fwrite(dzVol->fdata,sizeof(float),n,fpdz);
    fwrite(magVol->fdata,sizeof(float),n,fpmag);

    fclose(fpdx); fclose(fpdy); fclose(fpdz);  fclose(fpmag);
  } 

  return(1);
}    

int ReadSplatTable(ViewParameters *vpars,SplatTable *splat) {
 
  return(1);
}

void GetVolumeBoundingBox(Volume *vol,double isoVal)  {
  int nx,ny,nz,nxny;
  unsigned char *data;

  nx=vol->nx; ny=vol->ny; nz=vol->nz; nxny=nx*ny;  
  vol->bbLow.z=0; vol->bbHi.z=nz-1;
  vol->bbLow.y=0; vol->bbHi.y=ny-1;
  vol->bbLow.x=0; vol->bbHi.x=nx-1;
}

int ReadVolumeFile(char *volumeSetName,Volume *vol,double isoVal) {
  FILE *fp;
   char str[200];
   unsigned char *cdata;
  int i;
 
  if(!(fp=fopen(volumeSetName,"r")))
    return(0);
 
  
  fgets(str,200,fp);
  vol->nx=atoi(strtok(str," "));     
  vol->ny=atoi(strtok(NULL," "));     
  vol->nz=atoi(strtok(NULL," ")); 
    
  vol->n=vol->nx*vol->ny*vol->nz; 

  vol->fdata=(float*)calloc(vol->n,sizeof(float));  
   
  cdata=(unsigned char*)calloc(vol->n,sizeof(unsigned char));
  fread(cdata,sizeof(unsigned char),vol->n,fp);
  fclose(fp);
  for(i=0;i<vol->n;i++) { vol->fdata[i]=(float)cdata[i];}
     
  /*for(i=0;i<vol->n;i++){printf("\n %f",vol->fdata[i]);}   */
  GetVolumeBoundingBox(vol,isoVal);
  
  return(1);
}

void ReadSliceData(ViewParameters *vpars,Volume *vol){
}

void WriteImage(Image *img,int nxImg,int nyImg,int colImgFlag,char *outFile) {
  double max,min,scale,*inData;
  int x,y,n,i;
  unsigned char *imgData,r,g,b; 
  Color *colInData;
  FILE *fp;


  n=nxImg*nyImg;
  if(!colImgFlag)
  {
    inData=img->data;
    max=-10000.0;
    min=+10000.0;
    for(y=0;y<n;y++)
    {
      if(inData[y]>max)
        max=inData[y];
      if((inData[y]<min)&&(inData[y]>TINY))
        min=inData[y];    
    }

    imgData=(unsigned char*)calloc(nxImg*nyImg,sizeof(unsigned char));

    scale=255.0/(max-min);
    for(y=0;y<n;y++)
    {
      if(inData[y]<TINY)
        imgData[y]=0;
      else
        imgData[y]=(int)((inData[y]-min)*scale);
    }

    fp=fopen(outFile,"w");
    fprintf(fp,"P5\n%d %d\n255\n",nxImg,nyImg);
    fwrite(imgData,sizeof(unsigned char),nxImg*nyImg,fp);
    fclose(fp);
    free(imgData);
  }
  else  {
    fp=fopen(outFile,"w");
    fprintf(fp,"P5\n%d %d\n255\n",nxImg,nyImg);
    for(y=nyImg-1;y>=0;y--)
    {
      colInData=&img->colorData[y*nxImg];
      for(x=0;x<nxImg;x++)
      {
        r=colInData->r;
        g=colInData->g;;
        b=colInData->b;;
        fwrite(&r,sizeof(unsigned char),1,fp);
        colInData++;
      }
    }
    fclose(fp);
  }
}    

int CreateInterpolationKernel(InterpolationKernel *intKernel,
			      char *intKernelStr,double coeff,
			      double coeff2,ViewParameters *vpars) {
  int kernelType,nx,indx,ny,nz,ix,iy,iz,itmp,lim,numFiltWeights,j,i;
  double x,a1,b1,c1,d1,a2,b2,c2,d2,scale,scaleInv,*data,extent,u,v,w,nxny,t,
    tmp,smoothIntFilterCode,smoothIntFilterCoeff,
    smoothDerivFilterCode,smoothDerivFilterCoeff;
  float filtWeights[6];
  int sym=0;
  int center,index;
  smoothIntFilterCode   =vpars->smoothIntFilterCode;
  smoothIntFilterCoeff  =vpars->smoothIntFilterCoeff;
  smoothDerivFilterCode =vpars->smoothDerivFilterCode;
  smoothDerivFilterCoeff=vpars->smoothDerivFilterCoeff;
  if(!strcmp(intKernelStr,"smooth"))   {  
    sym=0;
    extent=3.0;
    nx=2*KERNEL_1D_NX*extent+1; ny=0; nz=0;
    data=(double*)calloc(nx,sizeof(double));
    kernelType=NON_SYMMETRIC;
    scale=(nx-1)/(2*extent);
    scaleInv=1./scale; 
    lim=(nx-1)/(2*extent);
    for(ix=0;ix<lim;ix++) 
    {  
      t=(float)ix/lim;
      numFiltWeights=getFiltWeights(smoothIntFilterCode,smoothIntFilterCoeff,
				    t,filtWeights);
      if(numFiltWeights==6)
      {
        for(j=0;j<6;j++)
          data[ix+j*lim]=filtWeights[j];
      }
      else
      {
        data[ix]=0;
        for(j=1;j<5;j++)
          data[ix+j*lim]=filtWeights[j-1];
        data[ix+5*lim]=0;
      }
    }        
    t=1.0;
    numFiltWeights=getFiltWeights(smoothIntFilterCode,smoothIntFilterCoeff,t
				  ,filtWeights);
    if(numFiltWeights==6)
      data[6*lim]=filtWeights[5];
    else if(numFiltWeights==5)
      data[4*lim]=filtWeights[4];
    /*   for(i=0;i<nx;i++)
      printf("%d %f\n",i,(float)data[i]);
    exit(1);*/
  }
  else if(!strcmp(intKernelStr,"smooth derivative"))
  {
    sym=0;
    extent=3.0;
    nx=2*KERNEL_1D_NX*extent+1; ny=0; nz=0;
    data=(double*)calloc(nx,sizeof(double));
    kernelType=SEPARABLE;
    scale=(nx-1.)/(2*extent);
    scaleInv=1./scale; 
    lim=(nx-1)/(2.*extent);
    for(ix=0;ix<lim;ix++) 
    {  
      t=(float)ix/lim;
      numFiltWeights=getFiltWeights(smoothDerivFilterCode,
				    smoothDerivFilterCoeff,t,filtWeights);
      if(numFiltWeights==6)
      {
        for(j=0;j<6;j++)
          data[ix+j*lim]=filtWeights[j];
      }
      else
      {
        data[ix]=0;
        for(j=1;j<5;j++)
          data[ix+j*lim]=filtWeights[j-1];
        data[ix+5*lim]=0;
      }
    }        
    t=1.0;
    numFiltWeights=getFiltWeights(smoothDerivFilterCode,
				  smoothDerivFilterCoeff,t,filtWeights);
    if(numFiltWeights==6)
      data[6*lim]=filtWeights[5];
    else if(numFiltWeights==5)
      data[4*lim]=filtWeights[4];
  }  

 /* for Mitchell-Netravali: coeff is -C, coeff2 is -B */ 

  else if((!strcmp(intKernelStr,"cubic spline"))||
	  (!strcmp(intKernelStr,"optimal")))
  {  
    sym=1;
    extent=2.0; 
    nx=KERNEL_1D_NX*extent+1; ny=0; nz=0;
    data=(double*)calloc(nx,sizeof(double));
    kernelType=SEPARABLE;
    scale=(nx-1)/extent;
    scaleInv=1./scale; 
    a1=coeff+2+(3.0/2.0)*coeff2;  a2=coeff+coeff2/6;
    b1=-(coeff+3+2*coeff2);       b2=-5*coeff-coeff2;
    c1=0;                         c2=8*coeff+2*coeff2;
    d1=1+(1.0/3.0)*coeff2;        d2=-4*coeff-(4.0/3.0)*coeff2;
    for(ix=0;ix<nx;ix++)
    {
      x=ix*scaleInv;
      if(x<1.0)
        data[ix]=(a1*x+b1)*x*x+d1;
      else
        data[ix]=((a2*x+b2)*x+c2)*x+d2;
    }
  }
  else if(!strcmp(intKernelStr,"cubic spline derivative"))
  {  
    sym=1;
    extent=2.0; 
    nx=KERNEL_1D_NX*extent+1; ny=0; nz=0;
    data=(double*)calloc(nx,sizeof(double));
    kernelType=SEPARABLE;
    scale=(nx-1)/extent;
    scaleInv=1./scale; 
    a1=6+3*coeff+(27.0/6.0)*coeff2;  a2=3*coeff+(1.0/2.0)*coeff2;
    b1=-(2*coeff+6+4*coeff2);        b2=-10*coeff-2*coeff2;
    c1=0;                            c2=8*coeff+2*coeff2;   
    for(ix=0;ix<nx;ix++)
    {
      x=ix*scaleInv;
      if(x<1.0)
        data[ix]=(a1*x+b1)*x;
      else
        data[ix]=(a2*x+b2)*x+c2;
    }
  }
  else if(!strcmp(intKernelStr,"trilinear"))
  {  
    sym=1;
    extent=1.0;
    nx=KERNEL_1D_NX*extent+1; ny=0; nz=0;
    scale=(nx-1)/extent;
    scaleInv=1.0/scale;
    data=(double*)calloc(nx,sizeof(double));
    kernelType=SEPARABLE;
    for(ix=0;ix<nx;ix++)
    {
      data[ix]=1.0-(ix*scaleInv);
    }
  }
  else if(!strcmp(intKernelStr,"optimal derivative"))
  {  
    sym=1;
    extent=2.0;
    nx=KERNEL_1D_NX*extent+1; ny=0; nz=0;
    scale=(nx-1.0)/extent;
    scaleInv=1.0/scale;
    data=(double*)calloc(nx,sizeof(double));
    kernelType=SEPARABLE;

    data[0]=0;
    data[(int)(KERNEL_1D_NX*extent/2.0)]=-0.5;
    data[(int)(KERNEL_1D_NX*extent)]=0;
    data[(int)(KERNEL_1D_NX*extent/2.0)]=-2./3.;
    data[(int)(KERNEL_1D_NX*extent)]=1./12.;
    itmp=KERNEL_1D_NX*extent/2.;
    for(ix=1;ix<itmp;ix++)
    {
      t=(ix*scaleInv);
      data[ix]=1.5*Sqr(t)-2.*t-0.5;
    }
    for(ix=itmp+1;ix<nx-1;ix++)
    {
      t=((ix-itmp)*scaleInv);
      data[ix]=-0.5*Sqr(t)+t-1./3.;
    }
    /*       for(ix=0;ix<nx;ix++)
      printf("\n%d %f",ix,data[ix]);*/
  }
  else 
    return(0);

  intKernel->extent=extent;
  intKernel->nx=nx; intKernel->ny=ny; intKernel->nz=nz;
  intKernel->scale=scale;
  intKernel->data=data;
  intKernel->type=kernelType;

   
  return(1);  
}

main(int argc, char *argv[]) {
  ViewParameters viewParams;
  Volume vol,dxVol,dyVol,dzVol,magVol,*volArr[6],shadeVol;
  SplatTable splatTable;
  clock_t time1;
  int i,j,n,filterCode;
  double *imageData;
  Color *colorImageData;
  char paramFileName[100];
  InterpolationKernel intKernel,derivIntKernel,derivKernel,*kernelArr[3];
  Image image;
   
  strcpy(paramFileName,"params");
  if(!ReadViewParameterFile(&viewParams,paramFileName)) {
    printf("\nCannot find view parameter file\n");
    exit(1);
  } 

  
  alphaDeriv=viewParams.derivKernelCoeffC;
  strcpy(normalEstimationMethod,viewParams.normalEstimationMethod);

  if(!strcmp(viewParams.splatMethod,"ray")||
     !strcmp(viewParams.splatMethod,"splat")) {
    if(!ReadSplatTable(&viewParams,&splatTable)) {
      printf("\nCannot find splat data file\n");
      exit(1);
    }
    viewParams.colImgFlag=FALSE;
  }
  else if(!strcmp(viewParams.splatMethod,"raycast")) {
    if(strcmp(viewParams.normalEstimationMethod,"FDH")&&
       strcmp(viewParams.normalEstimationMethod,"FHP")&&
       strcmp(viewParams.normalEstimationMethod,"SMO")) {
      printf("\nInterpolation menthod %s not known\n",
	     viewParams.normalEstimationMethod);
      exit(1);
    }

    if(!CreateInterpolationKernel(&intKernel,viewParams.intKernel,viewParams.
				  intKernelCoeffC,viewParams.intKernelCoeffB,
				  &viewParams)) {
      printf("\nCannot create interpolation kernel %s\n",viewParams.intKernel);
      exit(1);
    }
    strcpy(viewParams.derivIntKernel,viewParams.derivKernel);
    viewParams.derivIntKernelCoeffB=viewParams.derivKernelCoeffB;
    viewParams.derivIntKernelCoeffC=viewParams.derivKernelCoeffC;
    if(!CreateInterpolationKernel(&derivIntKernel,viewParams.derivIntKernel,
				  viewParams.derivIntKernelCoeffC,
				  viewParams.derivIntKernelCoeffB,
				  &viewParams)) {
      printf("\nCannot create interpolation kernel for gradient%s %s\n",
	     viewParams.derivIntKernel,viewParams.derivKernel);
      exit(1);
    }
    if(!strcmp(viewParams.normalEstimationMethod,"FHP")||
       !strcmp(viewParams.normalEstimationMethod,"SMO")) {
      sprintf(viewParams.derivKernel,"%s derivative",
	      viewParams.derivIntKernel);
      if(!CreateInterpolationKernel(&derivKernel,viewParams.derivKernel,
				    viewParams.derivKernelCoeffC,
				    viewParams.derivKernelCoeffB,
				    &viewParams)) {
        printf("\nCannot create derivative interpolation kernel %s\n",
	       viewParams.derivKernel);
        exit(1);
      }
    }
    viewParams.colImgFlag=TRUE;
    kernelArr[0]=&intKernel;
    kernelArr[1]=&derivKernel;
    kernelArr[2]=&derivIntKernel;
    int_filter_extent=kernelArr[0]->extent;
    der_filter_extent=kernelArr[2]->extent;

  }

  if(!strcmp(viewParams.volumeSetName,"nil")) {
    ReadSliceData(&viewParams,&vol);
  }
  else { 
    if(!ReadVolumeFile(viewParams.volumeSetName,&vol,viewParams.isovalue)) {
      printf("\nCannot find volume data datafile\n");
      exit(1);
    }
  }

  if(!strcmp(viewParams.splatMethod,"raycast")) {
    if(!strcmp(viewParams.normalEstimationMethod,"FDH")) {
      if(!CreateGradientVolumes(&vol,&dxVol,&dyVol,&dzVol,&magVol,&viewParams))
        exit(1);
    }
  }
  else if(!strcmp(viewParams.splatMethod,"ray")||
	  !strcmp(viewParams.splatMethod,"splat")) {
    if(!CreateGradientVolumes(&vol,&dxVol,&dyVol,&dzVol,&magVol,&viewParams))
      exit(1);
    shadeVol.nx=vol.nx; shadeVol.ny=vol.ny; shadeVol.nz=vol.nz; 
    shadeVol.data=(unsigned char*)calloc(vol.nx*vol.ny*vol.nz,
					 sizeof(unsigned char));
  }    
  
  volArr[0]=&vol; volArr[1]=&dxVol; volArr[2]=&dyVol; volArr[3]=&dzVol;
  volArr[4]=&magVol; volArr[5]=&shadeVol; 

  if(strstr(viewParams.splatMethod,"ray"))
  {
    viewParams.VRP.x+=(vol.nx-1)/2.0;   
    viewParams.VRP.y+=(vol.ny-1)/2.0;
    viewParams.VRP.z+=(vol.nz-1)/2.0;
    viewParams.lightSource.pos.x+=(vol.nx-1)/2.0;   
    viewParams.lightSource.pos.y+=(vol.ny-1)/2.0;
    viewParams.lightSource.pos.z+=(vol.nz-1)/2.0;
  }

  ComputeViewPlaneVectors(&viewParams);

  if(!viewParams.colImgFlag)
    imageData=(double*)calloc(viewParams.nxImg*
			      viewParams.nyImg,sizeof(double));
  else
    colorImageData=(Color*)calloc(viewParams.nxImg*viewParams.nyImg,
				  sizeof(Color));
  image.data=imageData;
  image.colorData=colorImageData;    

  

  clock();  
  for(i=0;i<viewParams.numberOfRuns;i++)
  {
    if(!strcmp(viewParams.splatMethod,"raycast"))
    {
      Raycast(&viewParams,volArr,colorImageData,kernelArr);
    }
  }

  time1=clock();
  strcpy(viewParams.outputFileName,"./Image.pgm");
  WriteImage(&image,viewParams.nxImg,viewParams.nyImg,
	     viewParams.colImgFlag,viewParams.outputFileName);
  printf("\n\nNeeded %9.8f s\n",time1/1000000.0);
 
  if(fork()==0) execlp("./XV",0);exit(0);
  /*Metrics(derivKernel.nx,derivKernel.data);*/
}