gridimplicit.cpp

隐式曲面求解的源代码

   double TM, VM, alpha;
   double discrim;

   v0=eval(v,x0,y0,z0); Assert(v0==v0);
   if(fabs(v0)<tol){
      if(!tHit) return true;
      TM=0;
      goto process_intersection;
   }
   v1=eval(v,x0+0.3333333333333*dt*ray.d[0],
             y0+0.3333333333333*dt*ray.d[1],
             z0+0.3333333333333*dt*ray.d[2]); Assert(v1==v1);
   if(fabs(v1)<tol){
      if(!tHit) return true;
      TM=0.3333333333333333;
      goto process_intersection;
   }
   v2=eval(v,x0+0.6666666666667*dt*ray.d[0],
             y0+0.6666666666667*dt*ray.d[1],
             z0+0.6666666666667*dt*ray.d[2]); Assert(v2==v2);
   if(fabs(v2)<tol){
      if(!tHit) return true;
      TM=0.6666666666666667;
      goto process_intersection;
   }
   v3=eval(v,x0+dt*ray.d[0],
             y0+dt*ray.d[1],
             z0+dt*ray.d[2]); Assert(v3==v3);
   if(fabs(v3)<tol){
      if(!tHit) return true;
      TM=1;
      goto process_intersection;
   }

   A=-4.5*v0+13.5*v1-13.5*v2+4.5*v3,
   B=9*v0-22.5*v1+18*v2-4.5*v3,
   C=-5.5*v0+9*v1-4.5*v2+v3,
   D=v0;  // coefficients of cubic: A*T^3+B*T^2+C*T+D  (where T=(t-t0)/(t1-t0))

#define SIGNCHANGE(a,b) (((a)>=0 && (b)<=0) || ((a)<=0 && (b)>=0))
   // look for first subinterval of [0,1] with a sign change in the cubic
   discrim=4*B*B-12*A*C;
   if(discrim<=0){ // if the cubic is strictly monotonic
      if(SIGNCHANGE(v0,v3)){
         T0=0; T1=1;
         V0=v0; V1=v3;
      }else
         return false;
   }else{ // we need to divide up the cubic into chunks to test for roots
      double rootDiscrim=sqrt(discrim), q;
      if(B<0) q=-.5*(B-rootDiscrim);
      else q=-.5*(B+rootDiscrim);
      double E1=q/A, E2=C/q;
      if(E2<E1) swap(E1,E2);
      // look for the intervals we need to check
      if(E1<0){
         if(E2<0){ // check [0,1]
            if(SIGNCHANGE(v0,v3)){
               T0=0; T1=1;
               V0=v0; V1=v3;
            }else
               return false;
         }else{
            if(E2<1){ // check [0,E2] and [E2,1]
               double E2v=D+E2*(C+E2*(B+E2*A));
               if(fabs(E2v)<tol){
                  if(!tHit) return true;
                  TM=E2;
                  goto process_intersection;
               }
               if(SIGNCHANGE(v0,E2v)){
                  T0=0; T1=E2;
                  V0=v0; V1=E2v;
               }else if(SIGNCHANGE(E2v,v3)){
                  T0=E2; T1=1;
                  V0=E2v; V1=v3;
               }else
                  return false;
            }else{ // check [0,1]
               if(SIGNCHANGE(v0,v3)){
                  T0=0; T1=1;
                  V0=v0; V1=v3;
               }else
                  return false;
            }
         }
      }else{
         if(E2<1){ // check [0,E1] [E1,E2] and [E2,1]
            double E1v=D+E1*(C+E1*(B+E1*A));
            if(fabs(E1v)<tol){
               if(!tHit) return true;
               TM=E1;
               goto process_intersection;
            }
            if(SIGNCHANGE(v0,E1v)){
               T0=0; T1=E1;
               V0=v0; V1=E1v;
            }else{
               double E2v=D+E2*(C+E2*(B+E2*A));
               if(fabs(E2v)<tol){
                  if(!tHit) return true;
                  TM=E2;
                  goto process_intersection;
               }
               if(SIGNCHANGE(E1v,E2v)){
                  T0=E1; T1=E2;
                  V0=E1v; V1=E2v;
               }else if(SIGNCHANGE(E2v,v3)){
                  T0=E2; T1=1;
                  V0=E2v; V1=v3;
               }else
                  return false;
            }
         }else{
            if(E1<1){ // check [0,E1] and [E1,1]
               double E1v=D+E1*(C+E1*(B+E1*A));
               if(fabs(E1v)<tol){
                  if(!tHit) return true;
                  TM=E1;
                  goto process_intersection;
               }
               if(SIGNCHANGE(v0,E1v)){
                  T0=0; T1=E1;
                  V0=v0; V1=E1v;
               }else if(SIGNCHANGE(E1v,v3)){
                  T0=E1; T1=1;
                  V0=E1v; V1=v3;
               }else
                  return false;
            }else{ // check [0,1]
               if(SIGNCHANGE(v0,v3)){
                  T0=0; T1=1;
                  V0=v0; V1=v3;
               }else
                  return false;
            }
         }
      }
   }

   // early exit if we don't care exactly where the intersection is
   if(!tHit) return true;

   // now do a few iterations of secant search
   if(V0<0){
      int i;
      for(i=0; i<5; ++i){
         alpha=V1/(V1-V0);
         TM=alpha*T0+(1-alpha)*T1;
         VM=D+TM*(C+TM*(B+TM*A));
         if(fabs(VM)<tol)
            break;
         else if(VM>0){
            T1=TM; V1=VM;
         }else{
            T0=TM; V0=VM;
         }
      }
   }else{
      int i;
      for(i=0; i<5; ++i){
         alpha=V1/(V1-V0);
         TM=alpha*T0+(1-alpha)*T1;
         VM=D+TM*(C+TM*(B+TM*A));
         if(fabs(VM)<tol)
            break;
         else if(VM<0){
            T1=TM; V1=VM;
         }else{
            T0=TM; V0=VM;
         }
      }
   }

   process_intersection:
   dt=dt*TM;
   *tHit=t0+dt;
   Assert(dg!=0);
   dg->p=ObjectToWorld(ray(*tHit));
   float xm=x0+dt*ray.d[0], ym=y0+dt*ray.d[1], zm=z0+dt*ray.d[2];
   float xr=1.f-xm, yr=1.f-ym, zr=1.f-zm;
   dg->nn[0]=yr*(zr*(v[4]-v[0])+zm*(v[5]-v[1]))
            +ym*(zr*(v[6]-v[2])+zm*(v[7]-v[3]));
   dg->nn[1]=xr*(zr*(v[2]-v[0])+zm*(v[3]-v[1]))
            +xm*(zr*(v[6]-v[4])+zm*(v[7]-v[5]));
   dg->nn[2]=xr*(yr*(v[1]-v[0])+ym*(v[3]-v[2]))
            +xm*(yr*(v[5]-v[4])+ym*(v[7]-v[6]));
   //dg->nn=Normal(ray(*tHit)-Point(0.5*dim[0],0.5*dim[1],0.5*dim[2]));
   dg->nn/=dg->nn.Length();
   dg->nn=ObjectToWorld(dg->nn);
   if(reverseOrientation)
      dg->nn*=-1.f;
   dg->u=0;
   dg->v=0;
   dg->shape=this;
   MakeTangentBasis(dg->nn, dg->dpdu, dg->dpdv);
   dg->dndu=Vector(0,0,0);
   dg->dndv=Vector(0,0,0);
   return true;
}

void GridImplicit::
GetShadingGeometry(const Transform &obj2world,
                   const DifferentialGeometry &dg,
                   DifferentialGeometry *dgShading) const
{
   *dgShading=dg;

   Point p=WorldToObject(dg.p);
   int i=Float2Int(p.x), j=Float2Int(p.y), k=Float2Int(p.z);
   // we rely on the user providing a band of empty cells around surface
   if(i<1) i=1; else if(i>dim[0]-3) i=dim[0]-3;
   if(j<1) j=1; else if(j>dim[1]-3) j=dim[1]-3;
   if(k<1) k=1; else if(k>dim[2]-3) k=dim[2]-3;
   // figure out coefficients for trilinear interpolation
   float xm=p.x-i, xr=1.f-xm;
   float ym=p.y-j, yr=1.f-ym;
   float zm=p.z-k, zr=1.f-zm;
   float c0=xr*yr*zr, c1=xr*yr*zm, c2=xr*ym*zr, c3=xr*ym*zm,
         c4=xm*yr*zr, c5=xm*yr*zm, c6=xm*ym*zr, c7=xm*ym*zm;
   dgShading->nn[0]=c0*(val(i+1,j,  k  ) - val(i-1,j,  k  ))
                   +c1*(val(i+1,j,  k+1) - val(i-1,j,  k+1))
                   +c2*(val(i+1,j+1,k  ) - val(i-1,j+1,k  ))
                   +c3*(val(i+1,j+1,k+1) - val(i-1,j+1,k+1))
                   +c4*(val(i+2,j,  k  ) - val(i,  j,  k  ))
                   +c5*(val(i+2,j,  k+1) - val(i,  j,  k+1))
                   +c6*(val(i+2,j+1,k  ) - val(i,  j+1,k  ))
                   +c7*(val(i+2,j+1,k+1) - val(i,  j+1,k+1));
   dgShading->nn[1]=c0*(val(i,  j+1,k  ) - val(i,  j-1,k  ))
                   +c1*(val(i,  j+1,k+1) - val(i,  j-1,k+1))
                   +c2*(val(i,  j+2,k  ) - val(i,  j,  k  ))
                   +c3*(val(i,  j+2,k+1) - val(i,  j,  k+1))
                   +c4*(val(i+1,j+1,k  ) - val(i+1,j-1,k  ))
                   +c5*(val(i+1,j+1,k+1) - val(i+1,j-1,k+1))
                   +c6*(val(i+1,j+2,k  ) - val(i+1,j,  k  ))
                   +c7*(val(i+1,j+2,k+1) - val(i+1,j,  k+1));
   dgShading->nn[2]=c0*(val(i,  j,  k+1) - val(i,  j,  k-1))
                   +c1*(val(i,  j,  k+2) - val(i,  j,  k  ))
                   +c2*(val(i,  j+1,k+1) - val(i,  j+1,k-1))
                   +c3*(val(i,  j+1,k+2) - val(i,  j+1,k  ))
                   +c4*(val(i+1,j,  k+1) - val(i+1,j,  k-1))
                   +c5*(val(i+1,j,  k+2) - val(i+1,j,  k  ))
                   +c6*(val(i+1,j+1,k+1) - val(i+1,j+1,k-1))
                   +c7*(val(i+1,j+1,k+2) - val(i+1,j+1,k  ));
   dgShading->nn/=dgShading->nn.Length();
   dgShading->nn=ObjectToWorld(dgShading->nn);
   if(reverseOrientation)
      dgShading->nn*=-1.f;
   MakeTangentBasis(dgShading->nn, dgShading->dpdu, dgShading->dpdv);
}

//============================== CreateShape =================================
extern "C" DLLEXPORT
Shape *CreateShape(const Transform &o2w,
                   bool reverseOrientation,
                   const ParamSet &params)
{
   string filename=params.FindOneString("filename","");
   return new GridImplicit(o2w, reverseOrientation, filename.c_str());
}