00001 #include <fstream.h>
00002 #include <math.h>
00003
00004 #include "marchingtetrahedra.h"
00005
00006
00007
00008
00009
00010
00011 vu1512121::vu1512121()
00012 {
00013 shading = 1;
00014 nTriangles = 0;
00015 printN = 1;
00016 }
00017
00018
00019
00020
00021
00022 vu1512121::vu1512121(const vu1512121& inst) : vu151212(inst)
00023 {
00024 }
00025
00026
00027
00028
00029
00030 vu1512121::~vu1512121()
00031 {
00032 }
00033
00034
00035
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
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
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
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
00085
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
00109
00110
00111 void vu1512121::render(void)
00112 {
00113 int p[4][3];
00114
00115
00116
00117 static const int index_offsets[12][4][3] = {{{0, 0, 1}, {1, 0, 1}, {1, 0, 0}, {1, 0, 2}},
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}},
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}},
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
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146 for (unsigned int x = 0; x < m_Dim1Size - 1; x++)
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++)
00151 {
00152 for (int i=0; i < 4; i++)
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);
00160 }
00161 }
00162 }
00163
00164
00165
00166
00167
00168 void vu1512121::DrawSurfaceInSimplex(int indices[4][3], float ht)
00169 {
00170
00171 float vt[4][3],nt[4][3] ;
00172 float g1[3], g2[3];
00173 int v1, v2;
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;
00188
00189 for (int m = 0; m < 4; m++)
00190 if (ht < values[m])
00191 whichSimplex |= (1 << m);
00192
00193 if ((whichSimplex == 0) || (whichSimplex == 15)) return;
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 {
00204 int edge = mtFacetData[whichSimplex][i];
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)
00231 {
00232 case WIRE:
00233 case FLAT:
00234 break;
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];
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;
00252 }
00253 break;
00254 }
00255 }
00256
00257
00258
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)
00302 {
00303 drawTriangle(vt[0], vt[2], vt[3], nt[0], nt[2], nt[3], (dot > 0.0));
00304 }
00305 }
00306
00307
00308
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++;
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);
00337
00338 glVertex3fv(v1);
00339 glVertex3fv(v2);
00340 glVertex3fv(v3);
00341 glEnd();
00342 }
00343 else
00344 {
00345 ComputeNormal(v1, v3, v2, n1);
00346
00347 glBegin(GL_TRIANGLES);
00348 glNormal3fv(n1);
00349
00350 glVertex3fv(v1);
00351 glVertex3fv(v3);
00352 glVertex3fv(v2);
00353 glEnd();
00354 }
00355 break;
00356 case SMOOTH:
00357 if (order)
00358 {
00359 glBegin(GL_TRIANGLES);
00360 glNormal3fv(n1);
00361 glVertex3fv(v1);
00362
00363 glNormal3fv(n2);
00364 glVertex3fv(v2);
00365
00366 glNormal3fv(n3);
00367 glVertex3fv(v3);
00368 glEnd();
00369 }
00370 else
00371 {
00372 glBegin(GL_TRIANGLES);
00373 glNormal3fv(n1);
00374 glVertex3fv(v1);
00375
00376 glNormal3fv(n3);
00377 glVertex3fv(v3);
00378
00379 glNormal3fv(n2);
00380 glVertex3fv(v2);
00381 glEnd();
00382 }
00383 break;
00384 }
00385 }
00386
00387
00388
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
00422
00423
00424 void vu1512121::ComputeNormal(float *p, float *q, float *r, float *n)
00425 {
00426 float pq[3], pr[3];
00427 float nNorm;
00428
00429 for (int i = 0; i < 3; i++)
00430 {
00431 pq[i] = q[i] - p[i];
00432 pr[i] = r[i] - p[i];
00433 }
00434
00435 n[0] = pq[1]*pr[2] - pr[1]*pq[2];
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];
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;
00448 }
00449 }
00450
00451
00452
00453
00454
00455 void vu1512121::InterpolatePoint(float x1, float y1, float z1, float h1,
00456 float x2, float y2, float z2, float h2,
00457 float ht, float* result)
00458 {
00459 float delta, one_delta;
00460
00461 if (h1 == h2)
00462 delta = 0.0;
00463 else
00464 delta = (h1 - ht) / (h1 - h2);
00465
00466 one_delta = 1.0 - delta;
00467
00468 result[0] = one_delta*x1 + delta*x2;
00469 result[1] = one_delta*y1 + delta*y2;
00470 result[2] = one_delta*z1 + delta*z2;
00471 }
00472
00473
00474
00475
00476
00477 void vu1512121::initOpenGL(void)
00478 {
00479 }