Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

marchingtetrahedra.cpp

Go to the documentation of this file.
00001 #include <fstream.h>
00002 #include <math.h>
00003 
00004 #include "marchingtetrahedra.h"
00005 
00006 
00007 //----------------------------------------------------------------------------
00008 //------------------------- The default constructor --------------------------
00009 //----------------------------------------------------------------------------
00010 
00011 vu1512121::vu1512121()
00012 {
00013   shading = 1;
00014   nTriangles = 0;
00015   printN = 1;
00016 }
00017 
00018 //----------------------------------------------------------------------------
00019 //------------------------- The copy constructor -----------------------------
00020 //----------------------------------------------------------------------------
00021 
00022 vu1512121::vu1512121(const vu1512121& inst) : vu151212(inst)
00023 {
00024 }
00025 
00026 //----------------------------------------------------------------------------
00027 //------------------------- The destructor -----------------------------------
00028 //----------------------------------------------------------------------------
00029 
00030 vu1512121::~vu1512121()
00031 {
00032 }
00033 
00034 //----------------------------------------------------------------------------
00035 //------------------------- The assignment operator --------------------------
00036 //----------------------------------------------------------------------------
00037 
00038 vu1512121& vu1512121::operator=(const vu1512121& rhs)
00039 {
00040     if (this != &rhs)
00041     {
00042         vu151212::operator=(rhs);
00043     }
00044     return *this;
00045 }
00046 
00047 //----------------------------------------------------------------------------
00048 //------------------------- public setViewVectors() --------------------------
00049 //----------------------------------------------------------------------------
00050 
00051 void vu1512121::setViewVectors(const vuVector& view,const vuVector& up,const vuVector& right)
00052 {
00053     m_View = view;
00054     m_Shift1 = right*3.2f;
00055     m_Shift2 = up*3.2f;
00056     m_Shift0 = (m_Shift1+m_Shift2)*(-0.5f);
00057 }
00058 
00059 
00060 //----------------------------------------------------------------------------
00061 //------------------------- public read() ------------------------------------
00062 //----------------------------------------------------------------------------
00063 
00064 bool vu1512121::read()
00065 {
00066     bool success = vu151212::read();
00067     if (!success) return false;
00068 
00069     return true;
00070 }
00071 
00072 //----------------------------------------------------------------------------
00073 //------------------------- public readRaw() ---------------------------------
00074 //----------------------------------------------------------------------------
00075 
00076 bool vu1512121::readRaw(void)
00077 {
00078     dword len;
00079     ifstream in;
00080 
00081 #ifdef IS_NOCREATE_NEEDED
00082     in.open(m_FileName, ios::in|ios::binary|ios::nocreate);
00083 #else
00084         // The nocreate is not available on the IRIX boxes that we were compiling
00085         // this code on, so we had to remove this part from
00086     in.open(m_FileName, ios::in|ios::binary);
00087 #endif
00088     if (!in.is_open()) return false;
00089 
00090     in.seekg(0, ios::end);
00091     len = in.tellg();
00092     in.seekg(0, ios::beg);
00093 
00094     in >> m_Dim1Size >> m_Dim2Size >> m_Dim3Size;
00095     if (in.fail()) return false;
00096     m_DataSize = m_Dim1Size*m_Dim2Size*m_Dim3Size;
00097 
00098     m_Data = new byte[m_DataSize];
00099     in.read((char *) (m_Data), int (m_DataSize));
00100     if (in.fail()) return false;
00101 
00102     in.close();
00103 
00104     return true;
00105 }
00106 
00107 //----------------------------------------------------------------------------
00108 //------------------------- public render() ----------------------------------
00109 //----------------------------------------------------------------------------
00110 
00111 void vu1512121::render(void)
00112 {
00113   int p[4][3];
00114   //  float rp[4][3];
00115 
00116   // bcc grid
00117   static const int index_offsets[12][4][3] = {{{0, 0, 1}, {1, 0, 1}, {1, 0, 0}, {1, 0, 2}}, // x-direction
00118                                                 {{0, 0, 1}, {1, 0, 1}, {1, 0, 2}, {1, 1, 2}},
00119                                                 {{0, 0, 1}, {1, 0, 1}, {1, 1, 2}, {1, 1, 0}},
00120                                                 {{0, 0, 1}, {1, 0, 1}, {1, 1, 0}, {1, 0, 0}},
00121                                                 {{0, 0, 1}, {0, 1, 1}, {0, 1, 0}, {0, 1, 2}}, // y-direction
00122                                                 {{0, 0, 1}, {0, 1, 1}, {0, 1, 2}, {1, 1, 2}},
00123                                                 {{0, 0, 1}, {0, 1, 1}, {1, 1, 2}, {1, 1, 0}},
00124                                                 {{0, 0, 1}, {0, 1, 1}, {1, 1, 0}, {0, 1, 0}},
00125                                                 {{0, 0, 1}, {0, 0, 3}, {0, 0, 2}, {0, 1, 2}}, // z-direction
00126                                                 {{0, 0, 1}, {0, 0, 3}, {0, 1, 2}, {1, 1, 2}},
00127                                                 {{0, 0, 1}, {0, 0, 3}, {1, 1, 2}, {1, 0, 2}},
00128                                                 {{0, 0, 1}, {0, 0, 3}, {1, 0, 2}, {0, 0, 2}}};
00129 
00130   /*
00131   // face divided
00132   static const int index_offsets[12][4][3] = {{{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
00133                                               {{1, 0, 0}, {0, 1, 0}, {1, 1, 0}, {0, 0, 1}},
00134                                               {{0, 0, 0}, {1, 0, 0}, {1, 0, 2}, {0, 0, 1}},
00135                                               {{0, 0, 0}, {0, 0, 2}, {1, 0, 2}, {0, 0, 1}},
00136                                               {{1, 0, 0}, {1, 1, 0}, {1, 1, 2}, {0, 0, 1}},
00137                                               {{1, 0, 0}, {1, 0, 2}, {1, 1, 2}, {0, 0, 1}},
00138                                               {{0, 1, 0}, {1, 1, 0}, {1, 1, 2}, {0, 0, 1}},
00139                                               {{0, 1, 0}, {0, 1, 2}, {1, 1, 2}, {0, 0, 1}},
00140                                               {{0, 0, 0}, {0, 1, 0}, {0, 1, 2}, {0, 0, 1}},
00141                                               {{0, 0, 0}, {0, 0, 2}, {0, 1, 2}, {0, 0, 1}},
00142                                               {{0, 0, 2}, {1, 0, 2}, {1, 1, 2}, {0, 0, 1}},
00143                                               {{0, 0, 2}, {0, 1, 2}, {1, 1, 2}, {0, 0, 1}}};
00144   */
00145 
00146   for (unsigned int x = 0; x < m_Dim1Size - 1; x++)           //        loop through cells
00147     for (unsigned int y = 0; y < m_Dim2Size - 1; y++)
00148         for (unsigned int z = 0; z < m_Dim3Size - 3; z+=2)
00149         {
00150           for (int s=0; s < 12; s++)  // 12 simplices
00151             {
00152               for (int i=0; i < 4; i++)  // calculate actual indices
00153                 {
00154                   p[i][0] = index_offsets[s][i][0] + x;
00155                   p[i][1] = index_offsets[s][i][1] + y;
00156                   p[i][2] = index_offsets[s][i][2] + z;
00157                 }
00158               
00159               DrawSurfaceInSimplex(p, threshold);                     //        call routine to draw the triangles
00160             }
00161         }
00162 }
00163 
00164 //----------------------------------------------------------------------------
00165 //------------------------- public DrawSurfaceInSimplex() ----------------------------------
00166 //----------------------------------------------------------------------------
00167 
00168 void vu1512121::DrawSurfaceInSimplex(int indices[4][3], float ht)       //      draws intersection with a given simplex
00169 {
00170   //SB: the changes I made in this file are just to prevent warnings...
00171   float vt[4][3],nt[4][3] /*, t[4][3] */;                               //      local array for the vertices
00172   float /*normal[3],*/ g1[3], g2[3];                    //      normal vector for triangles, gradients at vertices
00173   int v1, v2;                                           //      local copies of vertex ID
00174   float values[4];
00175 
00176   static const int mtFacetData[16][4] =
00177     {{0, 0, 0, 0}, {1, 3, 4, 0}, {1, 2, 5, 0}, {2, 3, 4, 5},
00178      {3, 2, 6, 0}, {1, 2, 6, 4}, {1, 5, 6, 3}, {4, 5, 6, 0},
00179      {4, 5, 6, 0}, {1, 3, 6, 5}, {1, 4, 6, 2}, {3, 6, 2, 0},
00180      {2, 5, 4, 3}, {1, 2, 5, 0}, {1, 3, 4, 0}, {0, 0, 0, 0}};
00181 
00182   static const int mtEdges[7][2] = {{0, 0}, {0, 1}, {1, 2}, {0, 2}, {0, 3}, {1, 3}, {2, 3}};
00183   
00184   for (int i = 0; i < 4; i++)
00185     values[i] = getDataValue(indices[i][0], indices[i][1], indices[i][2]);
00186  
00187   int whichSimplex = 0;                                 //      ID of type of cube for MC table
00188 
00189   for (int m = 0; m < 4; m++)                           //      loop through corners, computing facet lookup
00190     if (ht < values[m])                                 //      if the corner is above desired height
00191       whichSimplex |= (1 << m);                         //      set bit flag
00192 
00193   if ((whichSimplex == 0) || (whichSimplex == 15)) return;              //      bail out if no intersection
00194   
00195   int verts;
00196 
00197   if (mtFacetData[whichSimplex][3] == 0)
00198     verts = 3;
00199   else
00200     verts = 4;
00201 
00202   for (int i=0; i<verts; i++)
00203     { // for each vertex
00204       int edge = mtFacetData[whichSimplex][i];  //      local copies of edge ID
00205 
00206       v1 = mtEdges[edge][0];
00207       v2 = mtEdges[edge][1];
00208        
00209       float x1 = T*float(indices[v1][0]),
00210         y1 =     T*float(indices[v1][1]),
00211         z1 = 0.5*T*float(indices[v1][2]),
00212         x2 =     T*float(indices[v2][0]),
00213         y2 =     T*float(indices[v2][1]),
00214         z2 = 0.5*T*float(indices[v2][2]);
00215 
00216       if ((indices[v1][2] % 2) == 1)
00217         {
00218           x1 += 0.5*T;
00219           y1 += 0.5*T;
00220         }
00221 
00222       if ((indices[v2][2] % 2) == 1)
00223         {
00224           x2 += 0.5*T;
00225           y2 += 0.5*T;
00226         }
00227 
00228       InterpolatePoint(x1, y1, z1, values[v1], x2, y2, z2, values[v2], ht, vt[i]);
00229       
00230       switch(shading)  // compute normals
00231         {
00232         case WIRE:
00233         case FLAT:
00234           break;     // nothing to do
00235         case SMOOTH:
00236           computeGradient(indices[v1][0],
00237                           indices[v1][1],
00238                           indices[v1][2], g1, 1);
00239 
00240           computeGradient(indices[v2][0],
00241                           indices[v2][1],
00242                           indices[v2][2], g2, 1);
00243 
00244           InterpolatePoint(g1[0], g1[1], g1[2], values[v1], g2[0], g2[1], g2[2], values[v2], ht, nt[i]);
00245 
00246           double nNorm = nt[i][0]*nt[i][0] + nt[i][1]*nt[i][1] + nt[i][2]*nt[i][2];     //      compute the norm
00247       
00248           if (nNorm != 0.0)
00249             {
00250               nNorm = 1.0 / sqrt(nNorm);
00251               nt[i][0] *= nNorm; nt[i][1] *= nNorm; nt[i][2] *= nNorm;  //      and normalize the vector
00252             }
00253           break;
00254         }
00255     }
00256   
00257 
00258   // lets try to determine the order in one more way...
00259 
00260   int edge = mtFacetData[whichSimplex][0];
00261 
00262   v1 = mtEdges[edge][0];
00263   v2 = mtEdges[edge][1];
00264 
00265   float x1 = T*float(indices[v1][0]),
00266     y1 =     T*float(indices[v1][1]),
00267     z1 = 0.5*T*float(indices[v1][2]),
00268     x2 =     T*float(indices[v2][0]),
00269     y2 =     T*float(indices[v2][1]),
00270     z2 = 0.5*T*float(indices[v2][2]);
00271 
00272   if ((indices[v1][2] % 2) == 1)
00273     {
00274       x1 += 0.5*T;
00275       y1 += 0.5*T;
00276     }
00277   
00278   if ((indices[v2][2] % 2) == 1)
00279     {
00280       x2 += 0.5*T;
00281       y2 += 0.5*T;
00282     }
00283   
00284   float v[3], n[3];
00285 
00286   if (values[v1] < values[v2])
00287     {
00288       v[0] = x2 - x1; v[1] = y2 - y1; v[2] = z2 - z1;
00289     }
00290   else
00291     {
00292       v[0] = x1 - x2; v[1] = y1 - y2; v[2] = z1 - z2;
00293     }
00294 
00295   ComputeNormal(vt[0], vt[1], vt[2], n);
00296 
00297   float dot = v[0]*n[0] + v[1]*n[1] + v[2]*n[2];
00298 
00299   drawTriangle(vt[0], vt[1], vt[2], nt[0], nt[1], nt[2], (dot > 0.0));
00300 
00301   if (verts == 4)  // we have a quadrangle, so lets draw another triangle...
00302     {
00303       drawTriangle(vt[0], vt[2], vt[3], nt[0], nt[2], nt[3], (dot > 0.0));
00304     }
00305 } // DrawSimplex()
00306 
00307 //----------------------------------------------------------------------------
00308 //------------------------- public drawTriangle() ----------------------------------
00309 //----------------------------------------------------------------------------
00310 
00311 void vu1512121::drawTriangle(float v1[3], float v2[3], float v3[3],
00312                                float n1[3], float n2[3], float n3[3], bool order)
00313 {
00314   nTriangles++;        // count triangles
00315 
00316   switch(shading)
00317     {
00318     case WIRE:
00319       glBegin(GL_LINES);
00320       glVertex3fv(v1);
00321       glVertex3fv(v2);
00322           
00323       glVertex3fv(v2);
00324       glVertex3fv(v3);
00325           
00326       glVertex3fv(v3);
00327       glVertex3fv(v1);
00328       glEnd();
00329       break;
00330     case FLAT:
00331       if (order)
00332         {
00333           ComputeNormal(v1, v2, v3, n1);
00334           
00335           glBegin(GL_TRIANGLES);
00336           glNormal3fv(n1);                      //      set normal
00337           
00338           glVertex3fv(v1);                      //      set the first vertex
00339           glVertex3fv(v2);                      //      set the second vertex
00340           glVertex3fv(v3);                      //      set the third vertex
00341           glEnd();
00342         }
00343       else
00344         {
00345           ComputeNormal(v1, v3, v2, n1);
00346           
00347           glBegin(GL_TRIANGLES);
00348           glNormal3fv(n1);                      //      set normal
00349           
00350           glVertex3fv(v1);                      //      set the first vertex
00351           glVertex3fv(v3);                      //      set the second vertex
00352           glVertex3fv(v2);                      //      set the third vertex
00353           glEnd();
00354         }
00355       break;
00356     case SMOOTH:
00357       if (order)
00358         {
00359           glBegin(GL_TRIANGLES);
00360           glNormal3fv(n1);                      //      set normal
00361           glVertex3fv(v1);                      //      set the first vertex
00362           
00363           glNormal3fv(n2);                      //      and set it
00364           glVertex3fv(v2);                      //      set the second vertex
00365           
00366           glNormal3fv(n3);                      //      and set it
00367           glVertex3fv(v3);                      //      set the third vertex
00368           glEnd();
00369         }
00370       else
00371         {
00372           glBegin(GL_TRIANGLES);
00373           glNormal3fv(n1);                      //      set normal
00374           glVertex3fv(v1);                      //      set the first vertex
00375           
00376           glNormal3fv(n3);                      //      and set it
00377           glVertex3fv(v3);                      //      set the third vertex
00378           
00379           glNormal3fv(n2);                      //      and set it
00380           glVertex3fv(v2);                      //      set the second vertex
00381           glEnd();
00382         }
00383       break;
00384     }
00385 }
00386 
00387 //----------------------------------------------------------------------------
00388 //------------------------- public computeGradient() ------------------------------
00389 //----------------------------------------------------------------------------
00390 
00391 void vu1512121::computeGradient(int i, int j, int k, float *g, int mode)
00392 {
00393   switch(mode)
00394     {
00395     case 0:
00396       g[0] = -0.5*(getDataValue(i+1, j  , k  ) - getDataValue(i-1, j  , k  ));
00397       g[1] = -0.5*(getDataValue(i       , j+1, k  ) - getDataValue(i  , j-1, k  ));
00398       g[2] = -0.5*(getDataValue(i       , j  , k+2) - getDataValue(i  , j  , k-2));
00399       break;
00400     case 1:
00401       i += (k % 2);
00402       j += (k % 2);
00403 
00404       g[0] = ((getDataValue(i,j,k-1) - getDataValue(i-1, j, k-1)) +
00405         (getDataValue(i,j,k+1) - getDataValue(i-1, j, k+1)) +
00406         (getDataValue(i,j-1,k-1) - getDataValue(i-1, j-1, k-1)) +
00407         (getDataValue(i,j-1,k+1) - getDataValue(i-1, j-1, k+1)))/(4.0*T);
00408       g[1] = ((getDataValue(i,j,k-1) - getDataValue(i, j-1, k-1)) +
00409         (getDataValue(i,j,k+1) - getDataValue(i, j, k+1)) +
00410         (getDataValue(i-1,j,k-1) - getDataValue(i-1, j-1, k-1)) +
00411         (getDataValue(i-1,j,k+1) - getDataValue(i-1, j-1, k+1)))/(4.0*T);
00412       g[2] = ((getDataValue(i,j,k+1) - getDataValue(i, j, k-1)) +
00413         (getDataValue(i,j-1,k+1) - getDataValue(i, j-1, k-1)) +
00414         (getDataValue(i-1,j,k+1) - getDataValue(i-1, j, k-1)) +
00415         (getDataValue(i-1,j-1,k+1) - getDataValue(i-1, j-1, k-1)))/(4.0*T);
00416       break;
00417     }
00418 }
00419 
00420 //----------------------------------------------------------------------------
00421 //------------------------- public ComputeNormal() ------------------------------
00422 //----------------------------------------------------------------------------
00423 
00424 void vu1512121::ComputeNormal(float *p, float *q, float *r, float *n)           //      computes normal for triangle pqr, puts it in n
00425 {
00426   float pq[3], pr[3];                                                                                                                                   //      local variables for edge vectors
00427   float nNorm;                                                                                                                                          //      norm of this vector (I like unit vectors)
00428   
00429   for (int i = 0; i < 3; i++)                                                                                                                   //      loop to compute pq, pr
00430     {
00431       pq[i] = q[i] - p[i];
00432       pr[i] = r[i] - p[i];              //      just subtract . . .
00433     } // end of loop to compute vectors 
00434   
00435   n[0] = pq[1]*pr[2] - pr[1]*pq[2];                                                                                                     //      compute cross-product
00436   n[1] = pq[2]*pr[0] - pr[2]*pq[0];
00437   n[2] = pq[0]*pr[1] - pr[0]*pq[1];
00438 
00439   nNorm = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];                                                                                    //      compute the norm
00440 
00441   if (nNorm > 0.0)
00442     {
00443       nNorm = 1.0 / sqrt(nNorm);                                                                                                
00444       
00445       n[0] *= nNorm;
00446       n[1] *= nNorm;
00447       n[2] *= nNorm;                                                                            //      and adjust the vector
00448     }
00449 } // end of ComputeNormal()
00450 
00451 //----------------------------------------------------------------------------
00452 //------------------------- public InterpolatePoint() ------------------------------
00453 //----------------------------------------------------------------------------
00454 
00455 void vu1512121::InterpolatePoint(float x1, float y1, float z1, float h1,        //      linear interpolation along an edge
00456                                 float x2, float y2, float z2, float h2, 
00457                                 float ht, float* result)                                        
00458 { // InterpolatePoint()
00459   float delta, one_delta;                               //      parameters of interpolation
00460 
00461   if (h1 == h2)
00462     delta = 0.0;
00463   else
00464     delta = (h1 - ht) / (h1 - h2);              //      compute parameter from p1 to p2
00465 
00466   one_delta = 1.0 - delta;                                                                                                              //      and inverse parameter
00467 
00468   result[0] = one_delta*x1 + delta*x2;
00469   result[1] = one_delta*y1 + delta*y2;  //      compute resultant vertex
00470   result[2] = one_delta*z1 + delta*z2;                          
00471 } // InterpolatePoint() 
00472 
00473 //----------------------------------------------------------------------------
00474 //------------------------- public initOpenGL() ------------------------------
00475 //----------------------------------------------------------------------------
00476 
00477 void vu1512121::initOpenGL(void)
00478 {
00479 }

Generated on Wed Dec 15 21:20:29 2004 for vuVolume by  doxygen 1.3.9.1