#include "splat.h"

#define Bound(a) (a>1.0 ? 1.0 : a)
#define FHP 1
#define FDH 2
#define FHD 3
#define NIL 0
#define USE 1
#define CREATE 2

extern char normalEstimationMethod[100];
extern float alphaDeriv;

int ComputeVolumeIntersectionPoints(ViewParameters *vpars,Volume *vol,VolumeIntersectionPoints *volIntersectPts)
{
  Vector volVert[8],vtmp,vec1,vec2,normal,ray,vpu,vpv,imgOrg,nomu,nomv,
         faceVert[4],*ptr,denom,*ptrEntryV,*ptrExitV,*ptrExit,*ptrEntry;
  ImageSpaceVector projVert[4];
  int i,j,k,nxVol,nyVol,nzVol,npolys,nxImg,nyImg,n,uStart,uEnd,nxImg1,nyImg1,
      minVertIndx,vertIndxL,vertIndxR,nextVertIndxL,nextVertIndxR,vStart,u,v,
      *uMin,*uMax,vMin,vMax,vertIncrL,vertIncrR,uLim,
      bbLowx,bbLowy,bbLowz,bbHix,bbHiy,bbHiz;
  VolumeIntersectionPoints *volIntersectPtsPtr;
  unsigned char *imgData,ival;
  Polygon poly[6]=
       {{0,1,2,3},{0,4,5,1},{4,7,6,5},{7,3,2,6},{5,6,2,1},{0,3,7,4}};  
  double dot,nomuc,nomvc,dxImg,dyImg,nomMult,initValEntry,initValExit,t,
         min,mL,mR,vEnd,vEndL,vEndR,uL,uR,cux,cuy,cuz,cvx,cvy,cvz,cx,cy,cz,
         vx,vy,vz,ux,uy,uz,denominator,denomV,denomU,denomIncrv,denomIncru,
         EYEx,EYEy,EYEz,imgWidth1,imgHeight1;
  FILE *fp;
  
  vtmp.x=vol->nx/2.0-vpars->VRP.x;
  vtmp.y=vol->ny/2.0-vpars->VRP.y;
  vtmp.z=vol->nz/2.0-vpars->VRP.z;

  if(Dot(&vtmp,&vpars->VPN)<0)   /* looking away from the volume */
    return(0);    

  dxImg=vpars->imgWidth/vpars->nxImg;
  dyImg=vpars->imgHeight/vpars->nyImg;

  nxVol=vol->nx; nyVol=vol->ny; nzVol=vol->nz; 
  bbLowx=vol->bbLow.x; bbLowy=vol->bbLow.y; bbLowz=vol->bbLow.z;
  bbHix=vol->bbHi.x; bbHiy=vol->bbHi.y; bbHiz=vol->bbHi.z;
  
  nxImg=vpars->nxImg; nyImg=vpars->nyImg;
  
  if(vpars->stripIntegral)
  {
    imgWidth1=vpars->imgWidth;  
    imgHeight1=vpars->imgHeight;   
    nxImg1=nxImg+1;
    nyImg1=nyImg+1;
  }
  else
  {
    imgWidth1=vpars->imgWidth-dxImg;  
    imgHeight1=vpars->imgHeight-dyImg;   
    nxImg1=nxImg;
    nyImg1=nyImg;
  }
  
  imgData=(unsigned char *)calloc(nxImg1*nyImg1,1);

  npolys=6;
  
  volIntersectPts[0].pt=(Vector*)calloc(nxImg1*nyImg1,sizeof(Vector));
  volIntersectPts[1].pt=(Vector*)calloc(nxImg1*nyImg1,sizeof(Vector));

  uMin=(int*)malloc(nyImg1*sizeof(int));
  uMax=(int*)malloc(nyImg1*sizeof(int));

  for(i=0;i<nyImg1;i++)
  {
    uMin[i]=nxImg1;
    uMax[i]=-1;
  }
  vMin=nyImg1;
  vMax=-1;

  volVert[0].x=nxVol-1; volVert[0].y=0;        volVert[0].z=0;
  volVert[1].x=nxVol-1; volVert[1].y=0;        volVert[1].z=nzVol-1;
  volVert[2].x=nxVol-1; volVert[2].y=nyVol-1;  volVert[2].z=nzVol-1;
  volVert[3].x=nxVol-1; volVert[3].y=nyVol-1;  volVert[3].z=0;
  volVert[4].x=0;       volVert[4].y=0;        volVert[4].z=0;
  volVert[5].x=0;       volVert[5].y=0;        volVert[5].z=nzVol-1;
  volVert[6].x=0;       volVert[6].y=nyVol-1;  volVert[6].z=nzVol-1;
  volVert[7].x=0;       volVert[7].y=nyVol-1;  volVert[7].z=0;

/*  volVert[0].x=bbHix;  volVert[0].y=bbLowy; volVert[0].z=bbLowz;
  volVert[1].x=bbHix;  volVert[1].y=bbLowy; volVert[1].z=bbHiz;
  volVert[2].x=bbHix;  volVert[2].y=bbHiy;  volVert[2].z=bbHiz;
  volVert[3].x=bbHix;  volVert[3].y=bbHiy;  volVert[3].z=bbLowz;
  volVert[4].x=bbLowx; volVert[4].y=bbLowy; volVert[4].z=bbLowz;
  volVert[5].x=bbLowx; volVert[5].y=bbLowy; volVert[5].z=bbHiz;
  volVert[6].x=bbLowx; volVert[6].y=bbHiy;  volVert[6].z=bbHiz;
  volVert[7].x=bbLowx; volVert[7].y=bbHiy;  volVert[7].z=bbLowz;
*/
  imgOrg.x=vpars->VRP.x-vpars->VPU.x*(imgWidth1/2.0)-vpars->VPV.x*(imgHeight1/2.0);
  imgOrg.y=vpars->VRP.y-vpars->VPU.y*(imgWidth1/2.0)-vpars->VPV.y*(imgHeight1/2.0);
  imgOrg.z=vpars->VRP.z-vpars->VPU.z*(imgWidth1/2.0)-vpars->VPV.z*(imgHeight1/2.0);

  vpu.x=vpars->VPU.x*dxImg; vpu.y=vpars->VPU.y*dxImg; vpu.z=vpars->VPU.z*dxImg;
  vpv.x=vpars->VPV.x*dyImg; vpv.y=vpars->VPV.y*dyImg; vpv.z=vpars->VPV.z*dyImg;

  if(vpars->viewAngle==0)
  { 
    ray=vpars->VPN;

    nomu.x=vpv.y*ray.z-vpv.z*ray.y;
    nomu.y=vpv.z*ray.x-vpv.x*ray.z;  
    nomu.z=vpv.x*ray.y-vpv.y*ray.x;
    nomuc=imgOrg.x*(-nomu.x)+imgOrg.y*(-nomu.y)+imgOrg.z*(-nomu.z);

    nomv.x=vpu.z*ray.y-vpu.y*ray.z;
    nomv.y=vpu.x*ray.z-vpu.z*ray.x;  
    nomv.z=vpu.y*ray.x-vpu.x*ray.y;
    nomvc=imgOrg.x*(-nomv.x)+imgOrg.y*(-nomv.y)+imgOrg.z*(-nomv.z);

    denominator=vpv.x*nomv.x+vpv.y*nomv.y+vpv.z*nomv.z;

    for(i=0;i<npolys;i++)
    {
      faceVert[0]=volVert[poly[i].vertIndx[0]];
      faceVert[1]=volVert[poly[i].vertIndx[1]];
      faceVert[2]=volVert[poly[i].vertIndx[2]];
      faceVert[3]=volVert[poly[i].vertIndx[3]];

      vec1.x=faceVert[1].x-faceVert[0].x;
      vec1.y=faceVert[1].y-faceVert[0].y;
      vec1.z=faceVert[1].z-faceVert[0].z;
      vec2.x=faceVert[2].x-faceVert[1].x;
      vec2.y=faceVert[2].y-faceVert[1].y;
      vec2.z=faceVert[2].z-faceVert[1].z;

      normal=Cross(&vec2,&vec1);
      dot=Dot(&normal,&ray);
      if(dot<-TINY)  /* front face */
      {
        volIntersectPtsPtr=&volIntersectPts[0];
        vertIncrL=1;
        vertIncrR=3;
      }
      else if(dot>TINY)  /* back face */
      {
        volIntersectPtsPtr=&volIntersectPts[1];
        vertIncrL=3;
        vertIncrR=1;
      }
      else
        continue;
      for(j=0;j<4;j++)
      {
        projVert[j].u=(faceVert[j].x*nomu.x+faceVert[j].y*nomu.y+faceVert[j].z*nomu.z+nomuc)/denominator;
        projVert[j].v=(faceVert[j].x*nomv.x+faceVert[j].y*nomv.y+faceVert[j].z*nomv.z+nomvc)/denominator;
      }    
    
      if(fabs(normal.x)>TINY)
      {
        cux=0; cvx=0; cx=faceVert[0].x;
        cuy=vpu.y-vpu.x*ray.y/ray.x;
        cvy=vpv.y-vpv.x*ray.y/ray.x;
        cy=imgOrg.y+(faceVert[0].x-imgOrg.x)*ray.y/ray.x;
        cuz=vpu.z-vpu.x*ray.z/ray.x;
        cvz=vpv.z-vpv.x*ray.z/ray.x;
        cz=imgOrg.z+(faceVert[0].x-imgOrg.x)*ray.z/ray.x;
      }
      else if(fabs(normal.y)>TINY)
      {
        cuy=0; cvy=0; cy=faceVert[0].y;
        cux=vpu.x-vpu.y*ray.x/ray.y;
        cvx=vpv.x-vpv.y*ray.x/ray.y;
        cx=imgOrg.x+(faceVert[0].y-imgOrg.y)*ray.x/ray.y;
        cuz=vpu.z-vpu.y*ray.z/ray.y;
        cvz=vpv.z-vpv.y*ray.z/ray.y;
        cz=imgOrg.z+(faceVert[0].y-imgOrg.y)*ray.z/ray.y;
      }
      else 
      {
        cuz=0; cvz=0; cz=faceVert[0].z;
        cux=vpu.x-vpu.z*ray.x/ray.z;
        cvx=vpv.x-vpv.z*ray.x/ray.z;
        cx=imgOrg.x+(faceVert[0].z-imgOrg.z)*ray.x/ray.z;
        cuy=vpu.y-vpu.z*ray.y/ray.z;
        cvy=vpv.y-vpv.z*ray.y/ray.z;
        cy=imgOrg.y+(faceVert[0].z-imgOrg.z)*ray.y/ray.z;
      }

      min=GIANT;
      for(j=0;j<4;j++)
      {
        if(projVert[j].v<min)
        {
          min=projVert[j].v;
          minVertIndx=j;
        }      
      }
    
      if((int)projVert[(minVertIndx+vertIncrL)%4].v==(int)projVert[minVertIndx].v)
        vertIndxL=(minVertIndx+vertIncrL)%4; 
      else
        vertIndxL=minVertIndx;

      if((int)projVert[(minVertIndx+vertIncrR)%4].v==(int)projVert[minVertIndx].v)
        vertIndxR=(minVertIndx+vertIncrR)%4; 
      else
        vertIndxR=minVertIndx;
  
      nextVertIndxL=(vertIndxL+vertIncrL)%4;
      nextVertIndxR=(vertIndxR+vertIncrR)%4;
 
      mL=(projVert[nextVertIndxL].u-projVert[vertIndxL].u)/(projVert[nextVertIndxL].v-projVert[vertIndxL].v);
      mR=(projVert[nextVertIndxR].u-projVert[vertIndxR].u)/(projVert[nextVertIndxR].v-projVert[vertIndxR].v);

      vStart=Ceil(projVert[vertIndxL].v);
      if(vStart<0)
        vStart=0;

      vEndL=projVert[nextVertIndxL].v;
      vEndR=projVert[nextVertIndxR].v;

      vEnd=Min(vEndL,vEndR);
      if(vEnd>=nyImg1)
        vEnd=nyImg1-1;

/*      if(vEnd<vStart)
        continue; */

      uL=projVert[vertIndxL].u+mL*(vStart-projVert[vertIndxL].v);  
      uR=projVert[vertIndxR].u+mR*(vStart-projVert[vertIndxR].v);      

      vx=vStart*cvx+cx;
      vy=vStart*cvy+cy;
      vz=vStart*cvz+cz;
      while(1)
      {
        if(vStart<vMin)
         vMin=vStart;
        if(vEnd>vMax)
         vMax=vEnd;
        for(v=vStart;v<=vEnd;v++)
        {
          uStart=Ceil(uL);
          if(uStart<0)
            uStart=0;
          uEnd=(int)uR;
          if(uEnd>=nxImg1)
            uEnd=nxImg1-1;
          ptr=&(volIntersectPtsPtr->pt[v*nxImg1+uStart]);
          ux=vx+uStart*cux;
          uy=vy+uStart*cuy;
          uz=vz+uStart*cuz;
          if(uStart<uMin[v])
            uMin[v]=uStart;
          if(uEnd>uMax[v])
            uMax[v]=uEnd;
          for(u=uStart;u<=uEnd;u++)
          {
            ptr->x=ux; ptr->y=uy; ptr->z=uz;
            ux+=cux; uy+=cuy; uz+=cuz;
            ptr++;
          }
          uL+=mL;  uR+=mR;
          vx+=cvx;  vy+=cvy;  vz+=cvz;
        }
        if(vEnd==vEndL)
        {   
          vertIndxL=nextVertIndxL;
          nextVertIndxL=(vertIndxL+vertIncrL)%4;
          vEndL=projVert[nextVertIndxL].v;
          if(vEndL<=vEnd)
            break;
          vStart=Ceil(projVert[vertIndxL].v);
          if(vStart<0)
            vStart=0;
          mL=(projVert[nextVertIndxL].u-projVert[vertIndxL].u)/(projVert[nextVertIndxL].v-projVert[vertIndxL].v);
          uL=projVert[vertIndxL].u+mL*(vStart-projVert[vertIndxL].v);  
        }          
        else
        {   
          vertIndxR=nextVertIndxR;
          nextVertIndxR=(vertIndxR+vertIncrR)%4;
          vEndR=projVert[nextVertIndxR].v;
          if(vEndR<=vEnd)
            break;
          vStart=Ceil(projVert[vertIndxR].v);
          if(vStart<0)
            vStart=0;
          mR=(projVert[nextVertIndxR].u-projVert[vertIndxR].u)/(projVert[nextVertIndxR].v-projVert[vertIndxR].v);
          uR=projVert[vertIndxR].u+mR*(vStart-projVert[vertIndxR].v);  
        }          
        vEnd=Min(vEndL,vEndR);
        if(vEnd>=nyImg1)
          vEnd=nyImg1-1;
      }        
    }
  } 
  else
  {
    EYEx=vpars->EYE.x;  EYEy=vpars->EYE.y;  EYEz=vpars->EYE.z;
    for(i=0;i<8;i++)
    {
      volVert[i].x-=EYEx; volVert[i].y-=EYEy; volVert[i].z-=EYEz;
    }
    imgOrg.x-=EYEx;  imgOrg.y-=EYEy;  imgOrg.z-=EYEz;
    denom.x=vpu.z*vpv.y-vpu.y*vpv.z;
    denom.y=vpu.x*vpv.z-vpu.z*vpv.x;
    denom.z=vpu.y*vpv.x-vpu.x*vpv.y;

    nomu.x=vpv.z*imgOrg.y-vpv.y*imgOrg.z;
    nomu.y=vpv.x*imgOrg.z-vpv.z*imgOrg.x;
    nomu.z=vpv.y*imgOrg.x-vpv.x*imgOrg.y;

    nomv.x=vpu.y*imgOrg.z-vpu.z*imgOrg.y;
    nomv.y=vpu.z*imgOrg.x-vpu.x*imgOrg.z;
    nomv.z=vpu.x*imgOrg.y-vpu.y*imgOrg.x;

    for(i=0;i<npolys;i++)
    {
      faceVert[0]=volVert[poly[i].vertIndx[0]];
      faceVert[1]=volVert[poly[i].vertIndx[1]];
      faceVert[2]=volVert[poly[i].vertIndx[2]];
      faceVert[3]=volVert[poly[i].vertIndx[3]];

      vec1.x=faceVert[1].x-faceVert[0].x;
      vec1.y=faceVert[1].y-faceVert[0].y;
      vec1.z=faceVert[1].z-faceVert[0].z;
      vec2.x=faceVert[2].x-faceVert[1].x;
      vec2.y=faceVert[2].y-faceVert[1].y;
      vec2.z=faceVert[2].z-faceVert[1].z;

      ray=faceVert[1];

      normal=Cross(&vec2,&vec1);
      dot=Dot(&normal,&ray);
      if(dot<-TINY)  /* front face */
      {
        volIntersectPtsPtr=&volIntersectPts[0];
        vertIncrL=1;
        vertIncrR=3;
      }
      else if(dot>TINY)  /* back face */
      {
        volIntersectPtsPtr=&volIntersectPts[1];
        vertIncrL=3;
        vertIncrR=1;
      }
      else
        continue;
      for(j=0;j<4;j++)
      {
        denominator=(faceVert[j].x*denom.x+faceVert[j].y*denom.y+faceVert[j].z*denom.z);
        projVert[j].u=(faceVert[j].x*nomu.x+faceVert[j].y*nomu.y+faceVert[j].z*nomu.z)/denominator;
        projVert[j].v=(faceVert[j].x*nomv.x+faceVert[j].y*nomv.y+faceVert[j].z*nomv.z)/denominator;
      }    

      min=GIANT;
      for(j=0;j<4;j++)
      {
        if(projVert[j].v<min)
        {
          min=projVert[j].v;
          minVertIndx=j;
        }      
      }

      if((int)projVert[(minVertIndx+vertIncrL)%4].v==(int)projVert[minVertIndx].v)
        vertIndxL=(minVertIndx+vertIncrL)%4; 
      else
        vertIndxL=minVertIndx;

      if((int)projVert[(minVertIndx+vertIncrR)%4].v==(int)projVert[minVertIndx].v)
        vertIndxR=(minVertIndx+vertIncrR)%4; 
      else
        vertIndxR=minVertIndx;
  
      nextVertIndxL=(vertIndxL+vertIncrL)%4;
      nextVertIndxR=(vertIndxR+vertIncrR)%4;
 
      mL=(projVert[nextVertIndxL].u-projVert[vertIndxL].u)/(projVert[nextVertIndxL].v-projVert[vertIndxL].v);
      mR=(projVert[nextVertIndxR].u-projVert[vertIndxR].u)/(projVert[nextVertIndxR].v-projVert[vertIndxR].v);

      vStart=Ceil(projVert[vertIndxL].v);
      if(vStart<0)
        vStart=0;

      vEndL=projVert[nextVertIndxL].v;
      vEndR=projVert[nextVertIndxR].v;

      vEnd=Min(vEndL,vEndR);
      if(vEnd>=nyImg1)
        vEnd=nyImg1-1;

      uL=projVert[vertIndxL].u+mL*(vStart-projVert[vertIndxL].v);  
      uR=projVert[vertIndxR].u+mR*(vStart-projVert[vertIndxR].v);      

      if(fabs(normal.x)>TINY)
      {
        nomMult=faceVert[0].x;
        denomIncru=vpu.x;
        denomIncrv=vpv.x;
        denomV=imgOrg.x+vStart*denomIncrv;      
      }
      else if(fabs(normal.y)>TINY)
      {
        nomMult=faceVert[0].y;
        denomIncru=vpu.y;
        denomIncrv=vpv.y;
        denomV=imgOrg.y+vStart*denomIncrv;      
      }
      else 
      {
        nomMult=faceVert[0].z;
        denomIncru=vpu.z;
        denomIncrv=vpv.z;
        denomV=imgOrg.z+vStart*denomIncrv;      
      }

      cux=vpu.x*nomMult;
      cuy=vpu.y*nomMult;
      cuz=vpu.z*nomMult;
      cvx=vpv.x*nomMult;
      cvy=vpv.y*nomMult;
      cvz=vpv.z*nomMult;

      vx=imgOrg.x*nomMult+vStart*cvx;
      vy=imgOrg.y*nomMult+vStart*cvy;
      vz=imgOrg.z*nomMult+vStart*cvz;
      while(1)
      {
        if(vStart<vMin)
         vMin=vStart;
        if(vEnd>vMax)
         vMax=vEnd;
        for(v=vStart;v<=vEnd;v++)
        {
          uStart=Ceil(uL);
          if(uStart<0)
            uStart=0;
          uEnd=(int)uR;
          if(uEnd>=nxImg1)
            uEnd=nxImg1-1;
          denomU=denomV+uStart*denomIncru;
          ptr=&(volIntersectPtsPtr->pt[v*nxImg1+uStart]);
          ux=vx+uStart*cux;
          uy=vy+uStart*cuy;
          uz=vz+uStart*cuz;
          if(uStart<uMin[v])
            uMin[v]=uStart;
          if(uEnd>uMax[v])
            uMax[v]=uEnd;
          for(u=uStart;u<=uEnd;u++)
          {
            ptr->x=ux/denomU+EYEx; 
            ptr->y=uy/denomU+EYEy; 
            ptr->z=uz/denomU+EYEz;
/*            imgData[v*nxImg1+u]=ptr->x;*/
            ux+=cux; uy+=cuy; uz+=cuz;
            denomU+=denomIncru;
            ptr++;            
          }
          uL+=mL;  uR+=mR;
          vx+=cvx;  vy+=cvy;  vz+=cvz;
          denomV+=denomIncrv;
        }
        if(vEnd==vEndL)
        {   
          vertIndxL=nextVertIndxL;
          nextVertIndxL=(vertIndxL+vertIncrL)%4;
          vEndL=projVert[nextVertIndxL].v;
          if(vEndL<=vEnd)
            break;
          vStart=Ceil(projVert[vertIndxL].v);
          if(vStart<0)
            vStart=0;
          mL=(projVert[nextVertIndxL].u-projVert[vertIndxL].u)/(projVert[nextVertIndxL].v-projVert[vertIndxL].v);
          uL=projVert[vertIndxL].u+mL*(vStart-projVert[vertIndxL].v);  
        }          
        else
        {   
          vertIndxR=nextVertIndxR;
          nextVertIndxR=(vertIndxR+vertIncrR)%4;
          vEndR=projVert[nextVertIndxR].v;
          if(vEndR<=vEnd)
            break;
          vStart=Ceil(projVert[vertIndxR].v);
          if(vStart<0)
            vStart=0;
          mR=(projVert[nextVertIndxR].u-projVert[vertIndxR].u)/(projVert[nextVertIndxR].v-projVert[vertIndxR].v);
          uR=projVert[vertIndxR].u+mR*(vStart-projVert[vertIndxR].v);  
        }          
        vEnd=Min(vEndL,vEndR);
        if(vEnd>=nyImg1)
          vEnd=nyImg1-1;
      }               
    }
  }

  if(vpars->stripIntegral)
  {
    if(vpars->viewAngle==0)
    {
      vx=imgOrg.x; vy=imgOrg.y; vz=imgOrg.z;
      ray=vpars->VPN;
      ptrEntryV=volIntersectPts[0].pt;
      ptrExitV=volIntersectPts[1].pt;
      if(vpars->VPN.x>0)
      {
        initValEntry=nxVol+1;
        initValExit=-1;
      }
      else
      {
        initValEntry=-1;
        initValExit=nxVol+1;
      }
      for(v=0;v<nyImg1;v++)
      {
        uLim=uMin[v];
        ptrEntry=ptrEntryV;
        ptrExit=ptrExitV;
        ux=vx;  uy=vy;  uz=vz;
        for(u=0;u<uLim;u++)
        {
          ptrEntry->x=initValEntry;       
          t=(ptrEntry->x-ux)/ray.x;
          ptrEntry->y=uy+t*ray.y;
          ptrEntry->z=uz+t*ray.z;
          ptrExit->x=initValExit;         
          t=(ptrExit->x-ux)/ray.x;
          ptrExit->y=uy+t*ray.y;
          ptrExit->z=uz+t*ray.z;
          ux+=vpu.x;  uy+=vpu.y;  uz+=vpu.z;
          ptrEntry++;
          ptrExit++;
        }
        uLim=uMax[v]+1;
        ptrEntry=ptrEntryV+uLim;
        ptrExit=ptrExitV+uLim;
        ux=vx+uLim*vpu.x;  uy=vy+uLim*vpu.y;  uz=vz+uLim*vpu.z;
        for(u=uLim;u<nxImg1;u++)
        {
          ptrEntry->x=initValEntry;       
          t=(ptrEntry->x-ux)/ray.x;
          ptrEntry->y=uy+t*ray.y;
          ptrEntry->z=uz+t*ray.z;
          ptrExit->x=initValExit;         
          t=(ptrExit->x-ux)/ray.x;
          ptrExit->y=uy+t*ray.y;
          ptrExit->z=uz+t*ray.z;
          ux+=vpu.x;  uy+=vpu.y;  uz+=vpu.z;
          ptrEntry++;
          ptrExit++;
        }
        vx+=vpv.x;  vy+=vpv.y;  vz+=vpv.z;
        ptrEntryV+=nxImg1;
        ptrExitV+=nxImg1;
      }
    }
    else
    {
      vx=imgOrg.x; vy=imgOrg.y; vz=imgOrg.z;
      ptrEntryV=volIntersectPts[0].pt;
      ptrExitV=volIntersectPts[1].pt;
      if(vpars->VPN.x>0)
      {
        initValEntry=nxVol+1-EYEx;
        initValExit=-1-EYEx;
      }
      else
      {
        initValEntry=-1-EYEx;
        initValExit=nxVol+1-EYEx;
      }
      for(v=0;v<nyImg1;v++)
      {
        uLim=uMin[v];
        ptrEntry=ptrEntryV;
        ptrExit=ptrExitV;
        ux=vx;  uy=vy;  uz=vz;
        for(u=0;u<uLim;u++)
        {
          ptrEntry->x=initValEntry;       
          t=ptrEntry->x/ux;
          ptrEntry->x+=EYEx;
          ptrEntry->y=t*uy+EYEy;
          ptrEntry->z=t*uz+EYEz;
          ptrExit->x=initValExit;         
          t=ptrExit->x/ux;
          ptrExit->x+=EYEx;
          ptrExit->y=t*uy+EYEy;
          ptrExit->z=t*uz+EYEz;
          ux+=vpu.x;  uy+=vpu.y;  uz+=vpu.z;
          ptrEntry++;
          ptrExit++;
        }
        uLim=uMax[v]+1;
        ptrEntry=ptrEntryV+uLim;
        ptrExit=ptrExitV+uLim;
        ux=vx+uLim*vpu.x;  uy=vy+uLim*vpu.y;  uz=vz+uLim*vpu.z;
        for(u=uLim;u<nxImg1;u++)
        {
          ptrEntry->x=initValEntry;       
          t=ptrEntry->x/ux;
          ptrEntry->x+=EYEx;
          ptrEntry->y=t*uy+EYEy;
          ptrEntry->z=t*uz+EYEz;
          ptrExit->x=initValExit;         
          t=ptrExit->x/ux;
          ptrExit->x+=EYEx;
          ptrExit->y=t*uy+EYEy;
          ptrExit->z=t*uz+EYEz;
          ux+=vpu.x;  uy+=vpu.y;  uz+=vpu.z;
          ptrEntry++;
          ptrExit++;
        }
        vx+=vpv.x;  vy+=vpv.y;  vz+=vpv.z;
        ptrEntryV+=nxImg1;
        ptrExitV+=nxImg1;
      }
    }

/*    ptr=volIntersectPts[1].pt;
    for(v=0;v<nyImg1;v++)
      for(u=0;u<nxImg1;u++)
      {
        printf("\n%d\t%d\t%f\t%f\t%f",v,u,ptr->x,ptr->y,ptr->z);
        ptr++;
      }*/

    for(v=vMin;v<=vMax;v++)
    {
      if(uMin[v]>0)
        uMin[v]--;
      if(uMax[v]==nxImg1-1)
        uMax[v]--;  
    }
    if(vMin>0)
    {
      vMin--;
      uMin[vMin]=uMin[vMin+1];
      uMax[vMin]=uMax[vMin+1];
    }
    if(vMax==nyImg1-1)
      vMax--;
  }   

  volIntersectPts[0].vMax=vMax;
  volIntersectPts[0].vMin=vMin;
  volIntersectPts[0].uMax=uMax;
  volIntersectPts[0].uMin=uMin;
  volIntersectPts[0].nx=nxImg1;
  volIntersectPts[0].ny=nyImg1;

/*  for(i=0;i<nxImg1;i++)
    printf("\n%d %f",i,volIntersectPts[0].pt[50*nxImg+i].x);
*/
/*  fp=fopen("projImg.pgm","w");
  fprintf(fp,"P5\n%d %d\n255\n",nyImg1,nxImg1);
  n=nxImg1*nyImg1;
  for(i=0;i<n;i++)
  {
    ival=(int)volIntersectPts[0].pt[i].z+20;
    fwrite(&ival,1,1,fp);
  }
  fclose(fp);
*/
/*  fp=fopen("projImg.pgm","w");
  fprintf(fp,"P5\n%d %d\n255\n",nyImg1,nxImg1);
  fwrite(imgData,1,nxImg1*nyImg1,fp);
  fclose(fp);*/
  return(1);
}
  

void  Raycast(ViewParameters *vpars,Volume **volArr,Color *imageData,InterpolationKernel **kernelArr)
{  
  int nxImg,nxImg1,nyImg,nx,ny,nz,xVolIncr,yVolIncr,zVolIncr,uIncr,vIncr,
      nx1,ny1,nz1,vMax,vMin,*uMaxArr,*uMinArr,uMin,uMax,u,v,itmp,imgIndx,
      imgIndxV,xVolIndx,intKernelType,derivIntKernelType,nxIntKernel,
      nyIntKernel,nzIntKernel,nxnyIntKernel,ix,iy,iz,k,i,j,yprod,zprod,
      nxDerivIntKernel,nyDerivIntKernel,nzDerivIntKernel,nxnyDerivIntKernel,
      nsteps,xmin,xmax,ymin,ymax,zmin,zmax,volIndxX,volIndxZ,volIndxY,volIndx,
      normEstMethd,nImg,intImgFlag,firstHit,perspFlag,nxnyDerivKernel,
      nxDerivKernel,nyDerivKernel,nzDerivKernel,derivKernelType,indx,
      numAngle,cnt[3],
      totNumAnglesOverThresh,totNumAnglesUnderThresh,totNumEstAngles,
      numAnglesOverThresh,numAnglesUnderThresh,numEstAngles;
  Vector VPNabs,VPN,VRP,*entryPtsV,*exitPtsV,*entryPtsU,*exitPtsU,imgOrg,
         VPU,VPV,EYE,LGT,H,eye,lgt,normal,EYEP;
  Color lgtCol,objCol,ambientCol,colorSum;
  double tmp,rayx,rayy,rayz,absRayx,xEntry,xExit,sampleVal,convSumX,convSumZ,
         convSumY,dx,dy,dz,dydxRay,dzdxRay,imgPosx,imgPosxV,normMag,
         vpux,vpuy,vpuz,vpvx,vpvy,vpvz,intKernelExt,*intKernelData,
         intKernelScale,derivIntKernelExt,*derivIntKernelData,mag,
         derivIntKernelScale,isovalue,rayStepSize,tmp1,tmp2,tmp3,r,g,b,
         weight,totalWeight,totalWeightX,totalWeightZ,totalWeightY,x,y,z,dxImg,
         dyImg,ks,kd,ka,phongExp,dot,totalWeightx,totalWeighty,totalWeightz,
         convSumZx,convSumZy,convSumZz,convSumYx,convSumYy,convSumYz,
         totalWeightZx,totalWeightZy,totalWeightZz,totalWeightYx,totalWeightYy,
         totalWeightYz,sampleValx,sampleValy,sampleValz,*intImgX,*intImgY,
         *intImgZ,opacityLUT[255],opacitySum,opacity,gradMagLUT[255],
         derivKernelExt,*derivKernelData,derivKernelScale,dotprod,deltaLUT,
         expo,fx,fy,fz,xf,yf,zf,exactIsovalue,sampledIsovalue,angleError,
         totAngleError,angleThresh,f2x,f2y,f2z,f3x,f3y,f3z,fxsum,fysum,fzsum,
         fxEst,fyEst,fzEst,fxErrSum,fyErrSum,fzErrSum;
  VolumeIntersectionPoints volIntersectPts[2];
  FILE *fp,*fpIntImgX,*fpIntImgY,*fpIntImgZ,*fpErr;
  float *volData;
  float *dxData,*dyData,*dzData;
  Volume *vol,*dxVol,*dyVol,*dzVol;
  unsigned char *angleData;
  InterpolationKernel *derivKernel,*derivIntKernel,*intKernel;
  char str[100],filename[200];
  int cntfx,cntfy,cntfz,*histx,*histy,*histz;
  float taux,tauy,tauz,T,*localErrorx,*localErrory,*localErrorz,a1x,a1y,a1z;

  histx=(int*)calloc(501,sizeof(int));
  histy=(int*)calloc(501,sizeof(int));
  histz=(int*)calloc(501,sizeof(int));

  if(vpars->viewAngle>0)
    perspFlag=1;
  else
    perspFlag=0;

  vol=volArr[0];
  dxVol=volArr[1];
  dyVol=volArr[2];
  dzVol=volArr[3];

  intKernel=kernelArr[0];
  derivKernel=kernelArr[1];
  derivIntKernel=kernelArr[2];

  nxImg=vpars->nxImg;
  nxImg1=nxImg-1;
  nyImg=vpars->nyImg;
  nImg=nxImg*nyImg;

  angleData=(unsigned char*)calloc(nxImg*nyImg,1);

  xVolIncr=1;
  yVolIncr=vol->nx;
  zVolIncr=vol->ny*vol->nx;

  uIncr=1;
  vIncr=nxImg;

  VPNabs.x=fabs(vpars->VPN.x);
  VPNabs.y=fabs(vpars->VPN.y);
  VPNabs.z=fabs(vpars->VPN.z);
  if((VPNabs.y>VPNabs.x)||(VPNabs.z>VPNabs.x))
  {
    if(VPNabs.z>VPNabs.y)
    {
      Swap(vpars->VPN.z,vpars->VPN.x,tmp);
      Swap(vpars->VRP.z,vpars->VRP.x,tmp);
      Swap(vpars->VPU.z,vpars->VPU.x,tmp);
      Swap(vpars->VPV.z,vpars->VPV.x,tmp);
      Swap(vpars->lightSource.pos.z,vpars->lightSource.pos.x,tmp);
      Swap(vpars->EYE.z,vpars->EYE.x,tmp);
      Swap(xVolIncr,zVolIncr,itmp);
      Swap(vol->nz,vol->nx,itmp);
      Swap(vol->bbLow.z,vol->bbLow.x,itmp);
      Swap(vol->bbHi.z,vol->bbHi.x,itmp);
      vpars->VPV.x*=-1; vpars->VPV.y*=-1; vpars->VPV.z*=-1;
      vIncr*=-1;
    }
    else
    {
      Swap(vpars->VPN.y,vpars->VPN.x,tmp);
      Swap(vpars->VRP.y,vpars->VRP.x,tmp);
      Swap(vpars->VPU.y,vpars->VPU.x,tmp);
      Swap(vpars->VPV.y,vpars->VPV.x,tmp);
      Swap(vpars->lightSource.pos.y,vpars->lightSource.pos.x,tmp);
      Swap(vpars->EYE.y,vpars->EYE.x,tmp);
      Swap(vol->ny,vol->nx,itmp);
      Swap(vol->bbLow.y,vol->bbLow.x,itmp);
      Swap(vol->bbHi.y,vol->bbHi.x,itmp);
      Swap(xVolIncr,yVolIncr,itmp);
      vpars->VPU.x*=-1; vpars->VPU.y*=-1; vpars->VPU.z*=-1;
      uIncr*=-1;
    }
  }
  
  nx=vol->nx;
  ny=vol->ny;
  nz=vol->nz;
  nx1=nx-1;
  ny1=ny-1;
  nz1=nz-1;

  xVolIndx=xVolIncr;
  
  if(!ComputeVolumeIntersectionPoints(vpars,vol,volIntersectPts))
    return;

  VPN=vpars->VPN;
  VRP=vpars->VRP;
  VPU=vpars->VPU;
  VPV=vpars->VPV;

  EYE=VRP;
  LGT=vpars->lightSource.pos;
  lgtCol=vpars->lightSource.color;
  objCol=vpars->objProps.color;
  ambientCol=vpars->ambientColor;
  kd=vpars->objProps.kd;
  ks=vpars->objProps.ks;
  ka=vpars->objProps.ka;
  phongExp=vpars->objProps.phongExp;
  EYEP=vpars->EYE;

  dxImg=vpars->imgWidth/vpars->nxImg;
  dyImg=vpars->imgHeight/vpars->nyImg;

  vpux=VPU.x*dxImg;
  vpuy=VPU.y*dxImg;
  vpuz=VPU.z*dxImg;
  vpvx=VPV.x*dyImg;
  vpvy=VPV.y*dyImg;
  vpvz=VPV.z*dyImg;

  volData=vol->fdata;
  dxData=dxVol->fdata;
  dyData=dyVol->fdata;
  dzData=dzVol->fdata;

  if(!strcmp(vpars->normalEstimationMethod,"FHP"))
    normEstMethd=FHP;
  else if(!strcmp(vpars->normalEstimationMethod,"FDH"))
    normEstMethd=FDH;
  else if(!strcmp(vpars->normalEstimationMethod,"FHD"))
    normEstMethd=FHD;

  intKernelType=intKernel->type;  
  intKernelExt=intKernel->extent;
  nxIntKernel=intKernel->nx;
  nyIntKernel=intKernel->ny;
  nzIntKernel=intKernel->nz;
  nxnyIntKernel=nxIntKernel*nyIntKernel;
  intKernelData=intKernel->data;
  intKernelScale=intKernel->scale;

  derivIntKernelType=derivIntKernel->type;  
  derivIntKernelExt=derivIntKernel->extent;
  nxDerivIntKernel=derivIntKernel->nx;
  nyDerivIntKernel=derivIntKernel->ny;
  nzDerivIntKernel=derivIntKernel->nz;
  nxnyDerivIntKernel=nxDerivIntKernel*nyDerivIntKernel;
  derivIntKernelData=derivIntKernel->data;
  derivIntKernelScale=derivIntKernel->scale;

  if(normEstMethd==FHP)
  {
    derivKernelType=derivKernel->type;  
    derivKernelExt=derivKernel->extent;
    nxDerivKernel=derivKernel->nx;
    nyDerivKernel=derivKernel->ny;
    nzDerivKernel=derivKernel->nz;
    nxnyDerivKernel=nxDerivKernel*nyDerivKernel;
    derivKernelData=derivKernel->data;
    derivKernelScale=derivKernel->scale;
  }

  isovalue=vpars->isovalue;

  bzero(opacityLUT,255*sizeof(double));
  for(i=vpars->isovalue;i<255;i++)
    opacityLUT[i]=0.99+0.01*(i-vpars->isovalue)/(255.0-vpars->isovalue);

  for(i=0;i<255;i++)
    gradMagLUT[i]=1.0+i/100.0;

  if(!strcmp(vpars->intersectionImageFlag,"create"))
  {
    fpIntImgX=fopen("intImgX","w");
    intImgX=(double*)malloc(sizeof(double)*nImg);
    for(i=0;i<nImg;i++)
      intImgX[i]=-1; 
    intImgFlag=CREATE;
  }
  else if(!strcmp(vpars->intersectionImageFlag,"use"))
  {
    intImgFlag=USE;
    if((fpIntImgX=fopen("intImgX","r"))==NULL)
      intImgFlag=NIL;
    if(intImgFlag==NIL)
      printf("\nCould not find intersection images");
    else
    {
      intImgX=(double*)malloc(sizeof(double)*nImg);
      fread(intImgX,sizeof(double),nImg,fpIntImgX);  
      fclose(fpIntImgX);  
    }
  }
  else
    intImgFlag=NIL;

  rayx=VPN.x;
  rayy=VPN.y;
  rayz=VPN.z;
  dydxRay=rayy/rayx;
  dzdxRay=rayz/rayx;
  rayStepSize=vpars->rayStepSize;
  dx=rayx*rayStepSize;
  dy=rayy*rayStepSize;
  dz=rayz*rayStepSize;

  tmp1=(vpars->imgWidth-dxImg)/2.0;
  tmp2=(vpars->imgHeight-dyImg)/2.0;
  imgOrg.x=VRP.x-VPU.x*tmp1-VPV.x*tmp2;
  imgOrg.y=VRP.y-VPU.y*tmp1-VPV.y*tmp2;
  imgOrg.z=VRP.z-VPU.z*tmp1-VPV.z*tmp2;

  vMax=volIntersectPts[0].vMax;
  vMin=volIntersectPts[0].vMin;
  uMaxArr=volIntersectPts[0].uMax;
  uMinArr=volIntersectPts[0].uMin;

  entryPtsV=&volIntersectPts[0].pt[vMin*nxImg];
  exitPtsV=&volIntersectPts[1].pt[vMin*nxImg];
  
  imgPosxV=imgOrg.x+vMin*vpvx;
  if(vIncr<0)
    imgIndxV=(nyImg-1-vMin)*nxImg;
  else
    imgIndxV=vMin*nxImg;

  totAngleError=0;
  totNumEstAngles=0;
  totNumAnglesOverThresh=0;
  totNumAnglesUnderThresh=0;
  angleThresh=5;   /* for ml80 */
  angleThresh=25;  /* for ml40 */
  fxsum=0; fysum=0; fzsum=0;
  fxErrSum=0; fyErrSum=0; fzErrSum=0;
  cntfx=cntfy=cntfz=0;

  localErrorx=(float*)calloc(nxImg*nyImg,sizeof(float));
  localErrory=(float*)calloc(nxImg*nyImg,sizeof(float));
  localErrorz=(float*)calloc(nxImg*nyImg,sizeof(float));

  printf("\nraycast    "); 
  for(v=vMin;v<=vMax;v++)
  {
    sprintf(str,"\b\b\b%3d",v); fputs(str,stdout);
    uMax=uMaxArr[v];
    uMin=uMinArr[v];
    imgPosx=imgPosxV+vpux*uMin;
    entryPtsU=entryPtsV+uMin;
    exitPtsU=exitPtsV+uMin;
    if(uIncr<0)
      imgIndx=imgIndxV+(nxImg1-uMin);
    else
      imgIndx=imgIndxV+uMin;
    for(u=uMin;u<=uMax;u++)
    {    
      angleError=0;
      numEstAngles=0;
      numAnglesOverThresh=0;
      numAnglesUnderThresh=0;
      if(perspFlag)
      {
        rayx=entryPtsU->x-EYEP.x;
        rayy=entryPtsU->y-EYEP.y;
        rayz=entryPtsU->z-EYEP.z;
        mag=sqrt(Sqr(rayx)+Sqr(rayy)+Sqr(rayz));
        rayx/=mag; rayy/=mag; rayz/=mag; 

        dydxRay=rayy/rayx;
        dzdxRay=rayz/rayx;
        dx=rayx*rayStepSize;
        dy=rayy*rayStepSize;
        dz=rayz*rayStepSize;
      }

      if(intImgFlag!=USE)
        xEntry=Ceil(((entryPtsU->x-imgPosx)/dx))*dx+imgPosx;
      else
        xEntry=intImgX[imgIndx];

      xExit=Floor(((exitPtsU->x-imgPosx)/dx))*dx+imgPosx;
      nsteps=(xExit-xEntry)/dx+1;

      x=xEntry;
      y=entryPtsU->y+(xEntry-entryPtsU->x)*dydxRay;      
      z=entryPtsU->z+(xEntry-entryPtsU->x)*dzdxRay;
      /* printf("\n%d %d %f %f %f",v,u,x,y,z);      */
      colorSum.r=0; colorSum.g=0; colorSum.b=0; 
      opacitySum=0; 
      firstHit=1;
      for(k=0;k<nsteps;k++)
      {
        xmin=Max((Ceil(x-intKernelExt)),0);
        xmax=Min((Floor(x+intKernelExt)),nx1); 
        ymin=Max((Ceil(y-intKernelExt)),0);
        ymax=Min((Floor(y+intKernelExt)),ny1); 
        zmin=Max((Ceil(z-intKernelExt)),0);
        zmax=Min((Floor(z+intKernelExt)),nz1);
        volIndxZ=zmin*zVolIncr;
        sampleVal=0; 
        totalWeight=0;
        if(intKernelType==SEPARABLE)
	{         
          for(iz=zmin;iz<=zmax;iz++)
	  {
            volIndxY=volIndxZ+ymin*yVolIncr;
            convSumZ=0;
            totalWeightZ=0;
            for(iy=ymin;iy<=ymax;iy++)
            {
              volIndx=volIndxY+xmin*xVolIncr;
              convSumY=0; 
              totalWeightY=0;
              for(ix=xmin;ix<=xmax;ix++) 
	      {
                weight=intKernelData[(int)rint(fabs(x-ix)*intKernelScale)];
                totalWeightY+=weight;
                convSumY+=weight*volData[volIndx];
                volIndx+=xVolIncr;
  	      }
              weight=intKernelData[(int)rint(fabs(y-iy)*intKernelScale)];
              convSumZ+=convSumY*weight;
              totalWeightZ+=totalWeightY*weight;
              volIndxY+=yVolIncr;
            }
            weight=intKernelData[(int)rint(fabs(z-iz)*intKernelScale)];
            sampleVal+=convSumZ*weight;
            totalWeight+=totalWeightZ*weight;
            volIndxZ+=zVolIncr;
          }
        }
        else if(intKernelType==NON_SEPARABLE)
        {
          sampleVal=0; 
          for(iz=zmin;iz<=zmax;iz++)
          {
            volIndxY=volIndxZ+ymin*yVolIncr;
            zprod=(int)rint(fabs(z-iz)*intKernelScale)*nxnyIntKernel;
            for(iy=ymin;iy<=ymax;iy++)
            {
              yprod=(int)rint(fabs(y-iy)*intKernelScale)*nxIntKernel+zprod;
              volIndx=volIndxY+xmin*xVolIncr;
              for(ix=xmin;ix<=xmax;ix++) 
              {
                weight=intKernelData[yprod+(int)rint(fabs(x-ix)*intKernelScale)];
                sampleVal+=weight*volData[volIndx];
                totalWeight+=weight;
                volIndx++;
              }
              volIndxY+=yVolIncr;
            }
            volIndxZ+=zVolIncr;
          }          
        }
        if(sampleVal<0)
          sampleVal=0;
        sampleVal/=totalWeight;
        opacity=opacityLUT[(int)rint(sampleVal)];
	/*        if(sampleVal>isovalue)
          opacity=1.0;
        else
          opacity=0;*/        
        if(opacity>0)
	{           
          if((intImgFlag==CREATE)&&(firstHit))
	  {
            if(dx>0)
	    {
              intImgX[imgIndx]=x-Min(((int)(x/dx)),4)*dx;
            }
            else
	    {
              intImgX[imgIndx]=x-Min(((int)((nx1-x)/fabs(dx))),4)*dx;
            }   
            firstHit=0;
          }

          sampledIsovalue=sampleVal;
          r=ka*ambientCol.r*objCol.r;
          g=ka*ambientCol.g*objCol.g;
          b=ka*ambientCol.b*objCol.b;

          if(normEstMethd==FHP)
          { 
            xmin=Max((Ceil((x-derivIntKernelExt-0.001))),0);
            xmax=Min((Floor((x+derivIntKernelExt+0.001))),nx1); 
            ymin=Max((Ceil((y-derivIntKernelExt-0.001))),0);
            ymax=Min((Floor((y+derivIntKernelExt+0.001))),ny1); 
            zmin=Max((Ceil((z-derivIntKernelExt-0.001))),0);
            zmax=Min((Floor((z+derivIntKernelExt+0.001))),nz1);

            volIndxX=xmin*xVolIncr;
            sampleVal=0; 
            totalWeight=0;
            for(ix=xmin;ix<=xmax;ix++)
            {
              volIndxY=volIndxX+ymin*yVolIncr;
              convSumX=0;
              totalWeightX=0;
              for(iy=ymin;iy<=ymax;iy++)
              {
                volIndx=volIndxY+zmin*zVolIncr;
                convSumY=0; 
                totalWeightY=0;
                for(iz=zmin;iz<=zmax;iz++) 
                {
                  weight=derivIntKernelData[(int)rint(fabs(z-iz)*derivIntKernelScale)];
                  totalWeightY+=weight;
		  convSumY+=weight*volData[volIndx];
                  volIndx+=zVolIncr;
                }
                weight=derivIntKernelData[(int)rint(fabs(y-iy)*derivIntKernelScale)];
                convSumX+=convSumY*weight;
                totalWeightX+=totalWeightY*weight;
                volIndxY+=yVolIncr;
              }
              weight=derivKernelData[(int)rint(fabs(x-ix)*derivKernelScale)];
              if(x<ix) 
                weight*=-1;
              totalWeight+=totalWeightX*weight*(ix-x);
              sampleVal+=convSumX*weight;
              volIndxX+=xVolIncr;
            }
            normal.x=-sampleVal/totalWeight;
            volIndxY=ymin*yVolIncr;
            sampleVal=0; 
            totalWeight=0;
            for(iy=ymin;iy<=ymax;iy++)
	    {
              volIndxZ=volIndxY+zmin*zVolIncr;
              convSumY=0;
              totalWeightY=0;
              for(iz=zmin;iz<=zmax;iz++)
              {
                volIndx=volIndxZ+xmin*xVolIncr;
                convSumZ=0; 
                totalWeightZ=0;
                for(ix=xmin;ix<=xmax;ix++) 
                {
                  weight=derivIntKernelData[(int)rint(fabs(x-ix)*derivIntKernelScale)];
                  totalWeightZ+=weight;
                  convSumZ+=weight*volData[volIndx];
                  volIndx+=xVolIncr;
  	        }
                weight=derivIntKernelData[(int)rint(fabs(z-iz)*derivIntKernelScale)];
                convSumY+=convSumZ*weight;
                totalWeightY+=totalWeightZ*weight;
                volIndxZ+=zVolIncr;
              }
              weight=derivKernelData[(int)rint(fabs(y-iy)*derivKernelScale)];
              if(y<iy) 
                weight*=-1;
              totalWeight+=totalWeightY*weight*(iy-y);
              sampleVal+=convSumY*weight;
              volIndxY+=yVolIncr;
            }
            normal.y=-sampleVal/totalWeight;
            volIndxZ=zmin*zVolIncr;
            sampleVal=0; 
            totalWeight=0;
            for(iz=zmin;iz<=zmax;iz++)
            {
              volIndxY=volIndxZ+ymin*yVolIncr;
              convSumZ=0;
              totalWeightZ=0;
              for(iy=ymin;iy<=ymax;iy++)
              {
                volIndx=volIndxY+xmin*xVolIncr;
                convSumY=0; 
                totalWeightY=0;
                for(ix=xmin;ix<=xmax;ix++) 
                {
                  weight=derivIntKernelData[(int)rint(fabs(x-ix)*derivIntKernelScale)];
                  totalWeightY+=weight;
                  convSumY+=weight*volData[volIndx];
                  volIndx+=xVolIncr;
                }
                weight=derivIntKernelData[(int)rint(fabs(y-iy)*derivIntKernelScale)];
                convSumZ+=convSumY*weight; 
                totalWeightZ+=totalWeightY*weight;
                volIndxY+=yVolIncr;
              }
              weight=derivKernelData[(int)rint(fabs(z-iz)*derivKernelScale)];
              if(z<iz) 
                weight*=-1;
              totalWeight+=totalWeightZ*weight*(iz-z);
              sampleVal+=convSumZ*weight;
              volIndxZ+=zVolIncr;
            }
            normal.z=-sampleVal/totalWeight;
          }
          else if(normEstMethd==FDH)
          {
            xmin=Max((Ceil(x-derivIntKernelExt)),0);
            xmax=Min((Floor(x+derivIntKernelExt)),nx1); 
            ymin=Max((Ceil(y-derivIntKernelExt)),0);
            ymax=Min((Floor(y+derivIntKernelExt)),ny1); 
            zmin=Max((Ceil(z-derivIntKernelExt)),0);
            zmax=Min((Floor(z+derivIntKernelExt)),nz1);
            sampleValx=0; sampleValy=0; sampleValz=0; 
            totalWeightx=0; totalWeighty=0; totalWeightz=0;
            volIndxZ=zmin*zVolIncr;
            if(derivIntKernelType==SEPARABLE)
            {         
              for(iz=zmin;iz<=zmax;iz++)
              {
                volIndxY=volIndxZ+ymin*yVolIncr;
                convSumZx=0; convSumZy=0; convSumZz=0;
                totalWeightZx=0; totalWeightZy=0; totalWeightZz=0;
                for(iy=ymin;iy<=ymax;iy++)
                {
                  volIndx=volIndxY+xmin*xVolIncr;
                  convSumYx=0; convSumYy=0; convSumYz=0; 
                  totalWeightYx=0; totalWeightYy=0; totalWeightYz=0;
                  for(ix=xmin;ix<=xmax;ix++) 
                  {
                    weight=derivIntKernelData[(int)rint(fabs(x-ix)*derivIntKernelScale)];
                    totalWeightYx+=weight; 
                    totalWeightYy+=weight; 
                    totalWeightYz+=weight;
                    convSumYx+=weight*dxData[volIndx];
                    convSumYy+=weight*dyData[volIndx];
                    convSumYz+=weight*dzData[volIndx];
                    volIndx+=xVolIncr;
  	          }
                  weight=derivIntKernelData[(int)rint(fabs(y-iy)*derivIntKernelScale)];
                  convSumZx+=convSumYx*weight;
                  convSumZy+=convSumYy*weight;
                  convSumZz+=convSumYz*weight;
                  totalWeightZx+=totalWeightYx*weight;
                  totalWeightZy+=totalWeightYy*weight;
                  totalWeightZz+=totalWeightYz*weight;
                  volIndxY+=yVolIncr;
                }
                weight=derivIntKernelData[(int)rint(fabs(z-iz)*derivIntKernelScale)];
                sampleValx+=convSumZx*weight;
                sampleValy+=convSumZy*weight;
                sampleValz+=convSumZz*weight;
                totalWeightx+=totalWeightZx*weight;
                totalWeighty+=totalWeightZy*weight;
                totalWeightz+=totalWeightZz*weight;
                volIndxZ+=zVolIncr;
              }
            }
            else if(derivIntKernelType==NON_SEPARABLE)
            {
              for(iz=zmin;iz<=zmax;iz++)
              {
                volIndxY=volIndxZ+ymin*yVolIncr;
                zprod=(int)rint(fabs(z-iz)*derivIntKernelScale)*nxnyDerivIntKernel;
                for(iy=ymin;iy<=ymax;iy++)
                {    
                  yprod=(int)rint(fabs(y-iy)*intKernelScale)*nxDerivIntKernel+zprod;
                  volIndx=volIndxY+xmin*xVolIncr;
                  for(ix=xmin;ix<=xmax;ix++) 
                  {
                    weight=derivIntKernelData[yprod+(int)rint(fabs(x-ix)*derivIntKernelScale)];
                    sampleValx+=weight*dxData[volIndx];
                    sampleValy+=weight*dyData[volIndx];
                    sampleValz+=weight*dzData[volIndx];
                    totalWeightx+=weight;
                    totalWeighty+=weight;
                    totalWeightz+=weight;
                    volIndx++;
                  }
                  volIndxY+=yVolIncr;
                }
                volIndxZ+=zVolIncr;
              }          
            }
            normal.x=sampleValx/totalWeightx;  
            normal.y=sampleValy/totalWeighty;  
            normal.z=sampleValz/totalWeightz;
          }
	  /*         normal.x*=-1; normal.y*=-1; normal.z*=-1;*/

          normMag=sqrt(Sqr(normal.x)+Sqr(normal.y)+Sqr(normal.z));
          if(normMag<TINY)
            normMag=TINY; 
          fxEst=normal.x; fyEst=normal.y; fzEst=normal.z; 

          normal.x/=normMag; normal.y/=normMag; normal.z/=normMag;      

	  xf=(x-(nx-1)/2.0)*2.0/(nx-1);
          yf=(y-(ny-1)/2.0)*2.0/(ny-1);
          zf=(z-(nz-1)/2.0)*2.0/(nz-1);
          exactIsovalue=((1.-sin((M_PI*zf)/2.))+0.25*(1+cos(2*M_PI*6.*cos(M_PI*sqrt(xf*xf+yf*yf)/2.))))/2.5;
          fx=-(.6*sin(12.*M_PI*cos(0.5*sqrt(xf*xf+yf*yf)*M_PI))*M_PI*M_PI*sin(0.5*sqrt(xf*xf+yf*yf)*M_PI)*xf)/sqrt(xf*xf+yf*yf);
          fy=-(.6*sin(12.*M_PI*cos(0.5*sqrt(xf*xf+yf*yf)*M_PI))*M_PI*M_PI*sin(0.5*sqrt(xf*xf+yf*yf)*M_PI)*yf)/sqrt(xf*xf+yf*yf);
          fz=-(-.2*cos(.5*M_PI*zf)*M_PI);

          fxErrSum+=fabs(fxEst/(6.375)-fx); 
          fyErrSum+=fabs(fyEst/(6.375)-fy); 
          fzErrSum+=fabs(fzEst/(6.375)-fz);

          normMag=sqrt(Sqr(fx)+Sqr(fy)+Sqr(fz)); 
          fx/=normMag; fy/=normMag; fz/=normMag;

          dotprod=acos((fx*normal.x+fy*normal.y+fz*normal.z))*180.0/M_PI;
	  
	  /* derivatives 

          tmp1=0.5*M_PI*sqrt((xf*xf+yf*yf));
          tmp2=sin(12*M_PI*cos(tmp1));
          tmp3=cos(12.*M_PI*cos(tmp1));

          f2x=(.3*(-12.*M_PI*M_PI*cos(12.*M_PI*cos(tmp1))*xf*xf*sqrt(xf*xf+yf*yf)
               +sin(12.*M_PI*cos(tmp1))*M_PI*cos(tmp1)*xf*xf*sqrt(xf*xf+yf*yf)
               +12.*M_PI*M_PI*cos(12.*M_PI*cos(tmp1))*xf*xf*sqrt(xf*xf+yf*yf)*pow(cos(tmp1),2.)
               +2.*sin(12.*M_PI*cos(tmp1))*sin(tmp1)*yf*yf)*M_PI*M_PI)/pow(sqrt(xf*xf+yf*yf),1.5);

          f3x=(.15*(-144.*pow(xf*M_PI,4.)*tmp2*sin(tmp1)
               +144.*pow(xf*M_PI,4.)*tmp2*sin(tmp1)*pow(cos(tmp1),2.)
               -36.*pow(xf,4)*tmp3*pow(M_PI,3.)*sin(tmp1)*cos(tmp1)
               -pow(xf,4.)*tmp2*pow(M_PI,2.)*sin(tmp1)
               -pow(xf*M_PI*yf,2.)*tmp2*sin(tmp1)
               -36.*pow(xf*yf,2.)*tmp3*pow(M_PI,3.)*sin(tmp1)*cos(tmp1)
               -144.*pow(xf*yf,2.)*pow(M_PI,4.)*tmp2*sin(tmp1)
               +144.*pow(xf*yf*cos(tmp1),2.)*pow(M_PI,4.)*tmp2*sin(tmp1)
               +6.*tmp2*M_PI*cos(tmp1)*sqrt((xf*xf+yf*yf))*yf*yf  
               -72.*pow(M_PI*yf,2.)*tmp3*sqrt(xf*xf+yf*yf)
               +72.*pow(M_PI*yf*cos(tmp1),2.)*tmp3*sqrt(xf*xf+yf*yf)  
               -12.*tmp2*sin(tmp1)*yf*yf)*xf*M_PI*M_PI)/pow(xf*xf+yf*yf,5./2.);


          f2y=(.3*(-12.*M_PI*M_PI*cos(12.*M_PI*cos(tmp1))*yf*yf*sqrt(xf*xf+yf*yf)
               +sin(12.*M_PI*cos(tmp1))*M_PI*cos(tmp1)*yf*yf*sqrt(xf*xf+yf*yf)
               +12.*M_PI*M_PI*cos(12.*M_PI*cos(tmp1))*yf*yf*sqrt(xf*xf+yf*yf)*pow(cos(tmp1),2.)
               +2.*sin(12.*M_PI*cos(tmp1))*sin(tmp1)*xf*xf)*M_PI*M_PI)/pow(sqrt(xf*xf+yf*yf),1.5);

          f3y=(.15*(-144.*pow(yf*M_PI,4.)*tmp2*sin(tmp1)
               +144.*pow(yf*M_PI,4.)*tmp2*sin(tmp1)*pow(cos(tmp1),2.)
               -36.*pow(yf,4)*tmp3*pow(M_PI,3.)*sin(tmp1)*cos(tmp1)
               -pow(yf,4.)*tmp2*pow(M_PI,2.)*sin(tmp1)
               -pow(xf*M_PI*yf,2.)*tmp2*sin(tmp1)
               -36.*pow(xf*yf,2.)*tmp3*pow(M_PI,3.)*sin(tmp1)*cos(tmp1)
               -144.*pow(xf*yf,2.)*pow(M_PI,4.)*tmp2*sin(tmp1)
               +144.*pow(xf*yf*cos(tmp1),2.)*pow(M_PI,4.)*tmp2*sin(tmp1)
               +6.*tmp2*M_PI*cos(tmp1)*sqrt((xf*xf+yf*yf))*xf*xf  
               -72.*pow(M_PI*xf,2.)*tmp3*sqrt(xf*xf+yf*yf)
               +72.*pow(M_PI*xf*cos(tmp1),2.)*tmp3*sqrt(xf*xf+yf*yf)  
               -12.*tmp2*sin(tmp1)*xf*xf)*yf*M_PI*M_PI)/pow(xf*xf+yf*yf,5./2.);

          f2z=.1*sin(0.5*M_PI*zf)*M_PI*M_PI;
          f3z=.05*cos(0.5*M_PI*zf)*M_PI*M_PI*M_PI;

          tmp1=fabs(f3x/f2x);
          if(tmp1>500)
            histx[500]++;
          else
            histx[(int)(tmp1)]++;

          tmp1=fabs(f3y/f2y);
          if(tmp1>500)
            histy[500]++;
          else
            histy[(int)(tmp1)]++;

          tmp1=fabs(f3z/f2z);
          if(tmp1>500)
            histz[500]++;
          else
            histz[(int)(tmp1)]++;
            
          if(fabs(f2x)>TINY)
	  {
            fxsum+=fabs(f3x/f2x);
            cntfx++;
          }
          if(fabs(f2y)>TINY)
	  {
            fysum+=fabs(f3y/f2y);
            cntfy++;
          }
          if(fabs(f2z)>TINY)
	  {
            fzsum+=fabs(f3z/f2z);
            cntfz++;
          }

          tauz=z-(int)z;
          tauy=y-(int)y;
          taux=x-(int)x;
          T=2.0/40.;

          localErrorx[imgIndx]=fabs(fxEst/(6.375)-fx);
          localErrory[imgIndx]=fabs(fyEst/(6.375)-fy);
          localErrorz[imgIndx]=fabs(fzEst/(6.375)-fz);
	  */
	  /* 
          if((fabs(taux-0.5)>TINY)&&(fabs(taux)>TINY)&&(fabs(f2x)>TINY))
	  {
            localErrorx[imgIndx]=(f3x/f2x)*(T/12.)*((1+3*taux-3*taux*taux)/(taux*(1-taux)*(1-2*taux*taux)));
            if(localErrorx[imgIndx]>0)
              localErrorx[imgIndx]*=-1;
            localErrorx[imgIndx]-=0.5;
	  }       	      
          else
	  {
            tmp1=(3*taux-3*taux*taux)/(24*taux*taux*(taux-1)*(taux-1))-0.5;
            tmp2=-(2+3*taux-3*taux*taux)/(24*taux*taux*(taux-1)*(taux-1))-0.5;
            if(tmp2<tmp1)
              localErrorx[imgIndx]=tmp2;
            else
              localErrorx[imgIndx]=tmp1;
          }

          if((fabs(tauy-0.5)>TINY)&&(fabs(tauy)>TINY)&&(fabs(f2y)>TINY))
	  {
            localErrory[imgIndx]=(f3y/f2y)*(T/12.)*((1+3*tauy-3*tauy*tauy)/(tauy*(1-tauy)*(1-2*tauy*tauy)));
            if(localErrory[imgIndx]>0)
              localErrory[imgIndx]*=-1;
            localErrory[imgIndx]-=0.5;
	  }       	      
          else
	  {
            tmp1=(3*tauy-3*tauy*tauy)/(24*tauy*tauy*(tauy-1)*(tauy-1))-0.5;
            tmp2=-(2+3*tauy-3*tauy*tauy)/(24*tauy*tauy*(tauy-1)*(tauy-1))-0.5;
            if(tmp2<tmp1)
              localErrory[imgIndx]=tmp2;
            else
              localErrory[imgIndx]=tmp1;
          }

          if((fabs(tauz-0.5)>TINY)&&(fabs(tauz)>TINY)&&(fabs(f2z)>TINY))
	  {
            localErrorz[imgIndx]=(f3z/f2z)*(T/12.)*((1+3*tauz-3*tauz*tauz)/(tauz*(1-tauz)*(1-2*tauz*tauz)));
            if(localErrorz[imgIndx]>0)
              localErrorz[imgIndx]*=-1;
            localErrorz[imgIndx]-=0.5;
	  }       	      
          else
	  {
            tmp1=(3*tauz-3*tauz*tauz)/(24*tauz*tauz*(tauz-1)*(tauz-1))-0.5;
            tmp2=-(2+3*tauz-3*tauz*tauz)/(24*tauz*tauz*(tauz-1)*(tauz-1))-0.5;
            if(tmp2<tmp1)
              localErrorz[imgIndx]=tmp2;
            else
              localErrorz[imgIndx]=tmp1;
          }
	  */
                 
          /* for FDL 
        
          localErrorx[imgIndx]=fabs(f3x*T*T*T/6.*(1.+3*taux-3*taux*taux));
          localErrory[imgIndx]=fabs(f3y*T*T*T/6.*(1.+3*taux-3*taux*taux));
          localErrorz[imgIndx]=fabs(f3z*T*T*T/6.*(1.+3*taux-3*taux*taux));
	  */
          /* for FDH 

          localErrorx[imgIndx]=fabs(-f2x*T*T*(2.*alphaDeriv+1)*taux*(1-taux)*(1-2*taux));
          localErrory[imgIndx]=fabs(-f2y*T*T*(2.*alphaDeriv+1)*tauy*(1-tauy)*(1-2*tauy));
          localErrorz[imgIndx]=fabs(-f2z*T*T*(2.*alphaDeriv+1)*tauz*(1-tauz)*(1-2*tauz));
          if(localErrorx[imgIndx]<TINY)
            localErrorx[imgIndx]=fabs(f3x*T*T*T/6.*(1.+12.*(2.*alphaDeriv+1)*pow(taux*(taux-1),2.)));
          if(localErrory[imgIndx]<TINY)
            localErrory[imgIndx]=fabs(f3y*T*T*T/6.*(1.+12.*(2.*alphaDeriv+1)*pow(tauy*(tauy-1),2.)));
          if(localErrorz[imgIndx]<TINY)
            localErrorz[imgIndx]=fabs(f3z*T*T*T/6.*(1.+12.*(2.*alphaDeriv+1)*pow(tauz*(tauz-1),2.)));
	    */  
	  /* for FHP 
 
          a1x=T*(1.+(2.*alphaDeriv+1)*(-6*taux*taux+6*taux-1);
          a1y=T*(1.+(2.*alphaDeriv+1)*(-6*tauy*tauy+6*tauy-1);
          a1z=T*(1.+(2.*alphaDeriv+1)*(-6*tauz*tauz+6*tauz-1);

          localErrorx[imgIndx]=fabs(f2x*3*T*T*(2.*alphaDeriv+1)*taux*(1-taux)*(1-2*taux)/a1x);
          localErrory[imgIndx]=fabs(f2y*3*T*T*(2.*alphaDeriv+1)*tauy*(1-tauy)*(1-2*tauy)/a1y);
          localErrorz[imgIndx]=fabs(f2z*3*T*T*(2.*alphaDeriv+1)*tauz*(1-tauz)*(1-2*tauz)/a1z);
          if(localErrorx[imgIndx]<TINY)
            localErrorx[imgIndx]=fabs(f3x*T*T*T/6.*(1+12*alphaDeriv*taux-12*alphaDeriv*taux*taux-(2*alphaDeriv+1)*((((18*taux-36)*taux+24)*taux-6)*taux+1))/a1x);
          if(localErrory[imgIndx]<TINY)
            localErrory[imgIndx]=fabs(f3y*T*T*T/6.*(1+12*alphaDeriv*tauy-12*alphaDeriv*tauy*tauy-(2*alphaDeriv+1)*((((18*tauy-36)*tauy+24)*tauy-6)*tauy+1))/a1y);
          if(localErrorz[imgIndx]<TINY)
            localErrorz[imgIndx]=fabs(f3z*T*T*T/6.*(1+12*alphaDeriv*tauz-12*alphaDeriv*tauz*tauz-(2*alphaDeriv+1)*((((18*tauz-36)*tauz+24)*tauz-6)*tauz+1))/a1z);
	    */

	  /* printf("%f %f %f\n",localErrorx,localErrory,localErrorz);*/
          
	  /*          expo=6;
          xf=(x-nx/2.0)*10.0/nx;
          yf=(y-ny/2.0)*10.0/ny;
          zf=(z-nz/2.0)*10.0/nz;
          exactIsovalue=pow(fabs(xf),expo)+pow(fabs(yf),expo)+pow(fabs(zf),expo);
          fx=(fabs(xf)/xf)*expo*pow(fabs(xf),expo-1.0);
          fy=(fabs(yf)/yf)*expo*pow(fabs(yf),expo-1.0);
          fz=(fabs(zf)/zf)*expo*pow(fabs(zf),expo-1.0);
          normMag=sqrt(Sqr(fx)+Sqr(fy)+Sqr(fz)); 
          if(normMag<TINY)
            normMag=TINY; 
          fx/=normMag; fy/=normMag; fz/=normMag;
          dotprod=acos((fx*normal.x+fy*normal.y+fz*normal.z))*180.0/M_PI;
          if(v==50)
	  {
	    printf("\n%f %f %f %f %f %f",x,y,z,xf,yf,zf);
            printf("\n%d %f %f %f %f %f %f %f %f %f",u,sampledIsovalue,exactIsovalue,normal.x,fx,normal.y,fy,normal.z,fz,dotprod);
          }*/

          if(dotprod<0)
            dotprod*=-1;
 
          if(dotprod>angleThresh)
	  {
            numAnglesOverThresh++;
          }
          else
	  {
            numAnglesUnderThresh++;
          }
          angleError+=dotprod;
          numEstAngles++;

          lgt.x=LGT.x-x;  lgt.y=LGT.y-y;  lgt.z=LGT.z-z;
          NormalizeV(&lgt);
          if(kd>0)
          {
            dot=normal.x*lgt.x+normal.y*lgt.y+normal.z*lgt.z;          
            if(dot>0)
            {
              r+=objCol.r*dot*kd*lgtCol.r; 
              g+=objCol.g*dot*kd*lgtCol.g; 
              b+=objCol.b*dot*kd*lgtCol.b; 
            }
          }

          if(ks>0) 
          {
            eye.x=EYE.x-x;  eye.y=EYE.y-y;  eye.z=EYE.z-z;
            H.x=eye.x+lgt.x; H.y=eye.y+lgt.y; H.z=eye.z+lgt.z;
            NormalizeV(&H);
            dot=normal.x*H.x+normal.y*H.y+normal.z*H.z;
            if(dot>0)
            {
              dot=pow(dot,phongExp);          
              r+=dot*ks*lgtCol.r; g+=dot*ks*lgtCol.g; b+=dot*ks*lgtCol.b; 
            }
          }
	  /*opacity*=rayStepSize;*/
          if(r>1.0) r=1.0;
          if(g>1.0) g=1.0;
          if(b>1.0) b=1.0;
          colorSum.r=r*opacity*(1.0-opacitySum)+colorSum.r*opacitySum;
          colorSum.g=g*opacity*(1.0-opacitySum)+colorSum.g*opacitySum;
          colorSum.b=b*opacity*(1.0-opacitySum)+colorSum.b*opacitySum;
/*          if(normMag>255)
            printf("\n%f",normMag);
          opacity*=gradMagLUT[(int)rint(normMag)];*/
          opacitySum=opacity*(1.0-opacitySum)+opacitySum;
          if(opacitySum>0.999)
	  {
            imageData[imgIndx].r=255.0*Bound(colorSum.r);
            imageData[imgIndx].g=255.0*Bound(colorSum.g);
            imageData[imgIndx].b=255.0*Bound(colorSum.b);
            break;
          }
        }
        x+=dx; y+=dy; z+=dz;
      }
      if((intImgFlag==CREATE)&&(intImgX[imgIndx]==-1))
        intImgX[imgIndx]=x; 

      if(numAnglesOverThresh>0)
      {
        angleData[imgIndx]=255;
      }
      else
      {
        angleData[imgIndx]=(int)rint((angleError/numAnglesUnderThresh)*(255.0/angleThresh));
      }
      totNumAnglesOverThresh+=numAnglesOverThresh;
      totNumAnglesUnderThresh+=numAnglesUnderThresh;
      totNumEstAngles+=numEstAngles;
      totAngleError+=angleError;
 
      entryPtsU++;
      exitPtsU++;
      imgPosx+=vpux;
      imgIndx+=uIncr;
    }
    entryPtsV+=nxImg;
    exitPtsV+=nxImg;
    imgPosxV+=vpvx;
    imgIndxV+=vIncr;
  }

  /*  fpErr=fopen("errFile","w");
  fwrite(localErrorx,sizeof(float),nxImg*nyImg,fpErr);
  fwrite(localErrory,sizeof(float),nxImg*nyImg,fpErr);
  fwrite(localErrorz,sizeof(float),nxImg*nyImg,fpErr);
  fclose(fpErr);*/

  /*  printf("\nError: %f %5.3f %d %d\n",totAngleError/totNumEstAngles,100.0*totNumAnglesOverThresh/totNumEstAngles,totNumAnglesUnderThresh,totNumEstAngles);
  printf("\nDirectional errors: %f %f %f",fxErrSum/totNumEstAngles,fyErrSum/totNumEstAngles,fzErrSum/totNumEstAngles);*/
  /* printf("\nf3/f2: %f %f %f %d %d %d",fxsum/cntfx,fysum/cntfy,fzsum/cntfz,cntfx,cntfy,cntfz);*/

  /* for(i=0;i<=500;i++)
    printf("\n%d %d %d %d",i,histx[i],histy[i],histz[i]);*/

  /*printf("\n");
   sprintf(filename,"mlImages/ang_%s_%2.1f.pgm",normalEstimationMethod,alphaDeriv);
  printf("\n%s",filename);
  fp=fopen(filename,"w");*/

  fp=fopen("angle.pgm","w");
  fprintf(fp,"P5\n%d %d\n255\n",nxImg,nyImg);
  for(i=nyImg-1;i>=0;i--)
    fwrite(&angleData[i*nxImg],1,nxImg,fp);
  fclose(fp); 

  if(intImgFlag==CREATE)
  {
    fwrite(intImgX,sizeof(double),nImg,fpIntImgX);  
    fclose(fpIntImgX);  
  }
}