VolSlicer.cpp

Go to the documentation of this file.

//------------------------------------------------------------------------
//
//   Joe Kniss
//     6-17-03
//                   ________    ____   ___ 
//                  |        \  /    | /  /
//                  +---+     \/     |/  /
//                  +--+|  |\    /|     < 
//                  |  ||  | \  / |  |\  \ 
//                  |      |  \/  |  | \  \ 
//                   \_____|      |__|  \__\
//                       Copyright  2003 
//                      Joe Michael Kniss
//                   <<< jmk@cs.utah.edu >>>
//               "All Your Base are Belong to Us"
//-------------------------------------------------------------------------

//VolSlicer.cpp

/// A base class for slicing a volume.  More generaly, generating
///  volume samples.  A strategy for "VolRenBase". 
/// Also implements the standard volume slicer.

#include "VolSlicer.h"
#include <iostream>

using namespace gutz;
using namespace std;

///////////////////////////////////////////////////////////////////////////
/// 3D Slicing Vector Aligned
///////////////////////////////////////////////////////////////////////////

void VolSlicer::slice(mat4f        mv,   //model view matrix
                      VolytopeSP   vt,   //convex polytope with edge info
                      VolSamplesSP vs,   //storage for slice data
                      vec3f        axis) //axis to slice along eye-space coords
{

   arrayw1v3f vo = vt->getVerts();     //Verts of polytope in model-space
   arrayw1v3f tx = vt->getTcoords();   //Texture coordinates of polytope
   arrayo1v3f rv(vo.size(),vec3f(0)); //for the "Rotated Volume" (in eye-space)

   axis.normalize();  // can't be sure that the slice axis is a unit vector

   float maxval = -1000000000;  //used for testing for the
   float minval = +1000000000;  // nearest/farthes volume verticies

   //inverse model view matrix
   mat4f mvinv(mv.inv());       //inverse view matrix

   vec3f sn(mvinv.tdir(axis));  //slice plane normal in model-space
   sn.normalize();              // make sure it is a unit vector

   //reference point for slicing in eye-space
   vec3f eref(-axis*1000);

   //int  vnear = 0, vfar = 0;
   int  vnear,vfar;

   //translate model to world coords (view space) & find nearest/farthest points
   for(int i=0; i<vo.size(); ++i)
   {   
      rv[i] = mv.tpoint(vo[i]);           //Model -> eye space
      float dfr = (rv[i] - eref).norm();  // distance from reference
      if(maxval < dfr)                    // set if farthest point
      {
         maxval = dfr;
         vfar = i;
      }
      if(minval > dfr)                    // set if nearest point
      {                   
         minval = dfr;
         vnear = i; 
      }
   }

   if( (vnear == -1) || (vfar == -1) )
   {
     cerr << "VolSlicer::slice, invalid volume coordinates, check transforms that "
          << " generate the model view matrix!! " << endl;
     return;
   }

   //now compute the reference point in model-space
   vec3f mref = -sn * 1000;  //1000 out from (0,0,0) allong the negative slice plane normal

   //find the nearest slicing point (the closest sample)
   vec3f vec2near = (vo[vnear] - mref).dot(sn) * sn;    // the vector from reference to near along sn
   float sts = float((int)(vec2near.norm()/_sampleSpace)); // # of samples from reference to start
   vec3f nearPoint = sts*_sampleSpace * sn + mref;      // the sample nearest the reference point

   //find the farthest slicing point (same process as near point)
   vec3f vec2far = (vo[vfar] - mref).dot(sn) * sn;   // the vector from reference to far along sn
   float stf = float((int)(vec2far.norm()/_sampleSpace)+1);   // # of samples form reference to start
   vec3f farPoint = stf*_sampleSpace * sn + mref;         // the sample farthest from the reference point

   //distance to sample
   float dist = (nearPoint - farPoint).norm();

   //number of samples
   int samples = (int)(dist / _sampleSpace) + 1;

   //err(" #samples =", samples);

   //amount to move between samples
   vec3f delta = sn*_sampleSpace;

   //starting point for slicing
   vec3f sp = nearPoint;

   //slice data
   //vec3f poly[8];   // stores edge intersections (model-space)
   //vec3f tcoord[8]; // stores texture intersections (texture-space)

   arrayw1v2i edg = vt->getEdges();

   arrayo1v3f  veye(edg.size(),vec3f(0));  // stores transformed (eye-space) edge intersections

   arrayw1v3f   vout = vs->getVertAttrib()->getArray();
   arrayw1v4f   tout = vs->getTCoordAttrib()->getArray();
   arrayw1v3ui  iout = vs->getIdxAttrib()->getArray();
   arrayw1v2ui  pout = vs->getRangeAttrib()->getArray();

   int vStart = 0; //index of the first vertex in a new primitive
   int iStart = 0; //index of the first triangle index in a new primitive



   // could possibly benefit from optimization
   for(int i = 0 ; i < samples; ++i)
   { //for each slice

      sp += delta;

      int edges = 0; 
      //now check each edge of the volume for intersection with.. 
      //the plane defined by sp & sn

      //intersect the slice plane with the edges of our polytope
      for(int ne=0; ne < edg.size(); ++ne)
      {
         int v0 = edg[ne].v0;
         int v1 = edg[ne].v1;
         edges += intersect(vo[v0], vo[v1], tx[v0], tx[v1], rv[v0], rv[v1], sp, sn,
            vout[vStart+edges], tout[vStart+edges], veye[edges]);
      }

      //compute convex hull and sort, a little piece of magic from:
      // B.M.E. Moret & H.D. Shapiro "P to NP" pp. 453

      arrayo1i order(edges,0);
      //int *order = new int[edges]; //order of vertices
      vec3f ecen(0,0,0);            //center of all intersections (eye-space)
      vec3f mcen(0,0,0);            //           ''               (model-space)         
      vec3f tcen(0,0,0);            //           ''               (texture-space)

      for( int j=0; j<edges; j++)  //find the center of the points by averaging
      { 
         order[j] = j;           //set the initial order
         ecen += veye[j];       // & compute the centers
         mcen += vout[vStart+j];
         tcen += tout[vStart+j];
      }         
      ecen  *= 1.0f/(float)edges;
      mcen  *= 1.0f/(float)edges;
      tcen  *= 1.0f/(float)edges;

      for(int vi=0; vi<edges-1; vi++)  //for each vertex
      {          
         float theta = -10;            //find one with largest angle from center.. 
         int next = vi;               //the next nearest vertex
         for ( int k= vi; k<edges; k++)
         { 
            //... and check angle made between other edges
            vec3f dv(veye[order[k]] - ecen);
            if( (dv.x == ecen.x) && (dv.y == ecen.y))
            { //same as center? 
               next = k;
               //cout << "VolRenBase::Warning::degenerate slice" << endl;
               break; //out of this for-loop
            }

            //simply compute : y/(abs(x) + abs(y))
            float tt = (dv.y)/(mm_abs(dv.x) + mm_abs(dv.y));

            if( dv.x < 0.0f )      //check quadrants 2&3
            {
               tt = (float)(2.0f - tt);       
            }
            else if( dv.y < 0.0f ) //check quadrant 4
            {
               tt = (float)(4.0f + tt);       
            }

            if( theta < tt )
            {  //grab the max theta
               next = k;
               theta = tt;
            }
         } //end for(k) angle checking

         // I am not sure wich is worse: swapping 3 vectors for every edge
         // or: using an array to index into another array??? hmmm....
         //   should have payed more attention in class
         int tmp = order[vi];
         order[vi] = order[next];
         order[next] = tmp;               
      } //end for(j) edge /angle sort


      //now generate triangles for this slice, could be moved to the sort above
      //add the center of the slice to the verticies and texture coordinates
      vout[vStart+edges]   = mcen;
      tout[vStart+edges]   = tcen;

      //fill the index array based on the sort we did above
      for(int tri = 0; tri < edges; ++tri)
      {
         //iout[iStart + tri].x = vStart + order[tri];
         //iout[iStart + tri].y = vStart + order[((tri + 1) % edges)];
         //iout[iStart + tri].z = vStart + edges;

         iout[iStart + tri].z = vStart + order[tri];
         iout[iStart + tri].x = vStart + order[((tri + 1) % edges)];
         iout[iStart + tri].y = vStart + edges;
      }

      pout[i].v0 = iStart;
      pout[i].v1 = iStart + edges;

      vStart += edges + 1;
      iStart += edges;

   }// end for(i) each slice

   //err("num verts ", vStart);
   //err("num tris", iStart);

   vs->getVertAttrib()->setSize(vStart);
   vs->getTCoordAttrib()->setSize(vStart);
   vs->getIdxAttrib()->setSize(iStart);
   vs->getRangeAttrib()->setSize(samples);

}

//============================================================ INTERSECT
//======================================================================

int VolSlicer::intersect(
                         const vec3f &m0, const vec3f &m1, //line end points (model-space)
                         const vec3f &t0, const vec3f &t1, //texture for line end points (texture-space)
                         const vec3f &e0, const vec3f &e1, //view-space line end points (eye-space)
                         const vec3f &sp, const vec3f &sn, //plane point & plane norm  (in model-space)
                         vec3f &pnew,  vec4f &tnew, vec3f &vnew) //new values (respectively)
{
   float t = sn.dot(sp - m0) / sn.dot(m1 - m0);

   //note if the denominator is zero t is a NAN so we should have no problems?
   //probably not because NAN isn't between 0 & 1

   if( (t>=0) && (t<=1) )
   {
      pnew = m0 + t*(m1 - m0);
      tnew = t0 + t*(t1 - t0);
      vnew = e0 + t*(e1 - e0);
      return 1;
   }
   return 0;
} 

---

Send questions, comments, and bug reports to: