00001
00008 #ifndef SMOOTHMASK4D_H
00009 #define SMOOTHMASK4D_H
00010
00011 #ifdef __APPLE__
00012 #include <OpenGL/gl.h>
00013 #else
00014 #include <GL/gl.h>
00015 #endif
00016 #include <queue>
00017 #include <cmath>
00018 #include <error.h>
00019
00020 namespace Display
00021 {
00055 class SmoothMask4D
00056 {
00057 public:
00058 SmoothMask4D();
00059 ~SmoothMask4D();
00060
00066 template <typename DataType, typename MaskType>
00067 void update(const DataType* data,
00068 const MaskType* mask,
00069 unsigned int width,
00070 unsigned int height,
00071 unsigned int depth,
00072 unsigned int nb_frames);
00073
00075 void setFrame(unsigned int frame);
00076
00078 void display() const;
00079
00081 void clear();
00082
00083 private:
00084 template <typename DataType>
00085 void distanceTransform(const DataType* data,
00086 unsigned char* phi,
00087 std::vector<int> &band_index,
00088 unsigned int width,
00089 unsigned int height,
00090 unsigned int depth);
00091
00093 template <typename DataType>
00094 void march_cubes(const DataType* data,
00095 const unsigned char* mask,
00096 std::vector<int> &band_index,
00097 unsigned int width,
00098 unsigned int height,
00099 unsigned int depth);
00100
00103 struct XYZ { double x,y,z; };
00104 struct TRIANGLE { XYZ p[3]; };
00105 struct GRIDCELL { XYZ p[8]; double val[8]; };
00106 static XYZ VertexInterp(double isolevel, XYZ p1, XYZ p2, double valp1, double valp2);
00107 static int Polygonise(GRIDCELL grid,double isolevel,TRIANGLE *triangles);
00108
00109 private:
00110 GLuint m_meshid;
00111 unsigned int m_nb_frames;
00112 unsigned int m_frame;
00113 };
00114
00115 }
00116
00117
00118
00119 namespace Display
00120 {
00121
00122 template <typename DataType, typename MaskType>
00123 void SmoothMask4D::update(const DataType* data,
00124 const MaskType* mask,
00125 unsigned int width,
00126 unsigned int height,
00127 unsigned int depth,
00128 unsigned int nb_frames)
00129 {
00130
00131 if(glIsList(m_meshid)) glDeleteLists(m_meshid,m_nb_frames);
00132 m_nb_frames = nb_frames;
00133 m_meshid = glGenLists(m_nb_frames);
00134
00135
00136 unsigned char* phi = new unsigned char [width*height*depth];
00137
00138 for(unsigned int frame = 0; frame<m_nb_frames; ++frame)
00139 {
00140 glNewList(m_meshid+frame, GL_COMPILE);
00141
00142
00143 glMatrixMode(GL_MODELVIEW);
00144 glPushMatrix();
00145 glTranslatef(-0.5f,-0.5f,-0.5f);
00146 glScalef(1.0f/width, 1.0f/height, 1.0f/depth);
00147
00148
00149 glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
00150 glEnable(GL_COLOR_MATERIAL);
00151
00152
00153 glEnable(GL_LIGHTING);
00154 glEnable(GL_LIGHT0);
00155
00156
00157 std::vector<int> band_index;
00158 distanceTransform(&mask[frame*width*height*depth], phi, band_index, width, height, depth);
00159
00160
00161 if(data)
00162 march_cubes(&data[frame*width*height*depth], phi, band_index, width, height, depth);
00163 else
00164 march_cubes(&mask[frame*width*height*depth], phi, band_index, width, height, depth);
00165
00166
00167 glDisable(GL_LIGHT0);
00168 glDisable(GL_LIGHTING);
00169 glDisable(GL_COLOR_MATERIAL);
00170
00171
00172 glPopMatrix();
00173 glEndList();
00174 }
00175
00176
00177 delete [] phi;
00178 }
00179
00180
00181 template <typename DataType>
00182 void SmoothMask4D::distanceTransform(const DataType* mask,
00183 unsigned char* phi,
00184 std::vector<int> &band_index,
00185 unsigned int width,
00186 unsigned int height,
00187 unsigned int depth)
00188 {
00189 int neighbors[][3] = {
00190 {0,0,-1}, {0,-1,0}, {-1,0,0},
00191 {0,0,+1}, {0,+1,0}, {+1,0,0}
00192 };
00193
00195
00196
00197
00198 std::queue<int> inner_boundary[2];
00199 std::queue<int> outer_boundary[2];
00200 bool current_list = 0;
00201 unsigned char zero = 128;
00202 const int max_radius = 5;
00203
00204
00205 for(unsigned int i=0; i<width*height*depth; ++i) phi[i] = zero;
00206
00207
00208 for(unsigned int k=0; k<depth; ++k)
00209 for(unsigned int j=0; j<height; ++j)
00210 for(unsigned int i=0; i<width; ++i)
00211 {
00212 int index = (k*height+j)*width+i;
00213 if(mask[index] != 0)
00214 {
00215 for(unsigned int m=0; m<6; ++m)
00216 {
00217 int ii = i + neighbors[m][0];
00218 int jj = j + neighbors[m][1];
00219 int kk = k + neighbors[m][2];
00220
00221 if(ii >=0 && ii < (int)width && jj >= 0 && jj < (int)height && kk>=0 && kk<(int)depth)
00222 {
00223
00224 if(mask[(kk*height+jj)*width+ii] == 0)
00225 {
00226 phi[index] = zero+1;
00227 inner_boundary[current_list].push(index);
00228 break;
00229 }
00230 }
00231 }
00232 }
00233 }
00234
00235
00236 for(unsigned int k=0; k<depth; ++k)
00237 for(unsigned int j=0; j<height; ++j)
00238 for(unsigned int i=0; i<width; ++i)
00239 {
00240 int index = (k*height+j)*width+i;
00241 if(mask[index] == 0)
00242 {
00243 for(unsigned int m=0; m<6; ++m)
00244 {
00245 int ii = i + neighbors[m][0];
00246 int jj = j + neighbors[m][1];
00247 int kk = k + neighbors[m][2];
00248
00249 if(ii >=0 && ii < (int)width && jj >= 0 && jj < (int)height && kk>=0 && kk<(int)depth)
00250 {
00251
00252 if(mask[(kk*height+jj)*width+ii] != 0)
00253 {
00254 phi[index] = zero-1;
00255 outer_boundary[current_list].push(index);
00256 break;
00257 }
00258 }
00259 }
00260 }
00261 }
00262
00263
00264 current_list = 0;
00265 int layer = 1;
00266 while(!inner_boundary[current_list].empty() && layer < max_radius)
00267 {
00268
00269 while(!inner_boundary[current_list].empty())
00270 {
00271
00272 int index = inner_boundary[current_list].front();
00273 inner_boundary[current_list].pop();
00274
00275
00276 if(phi[index] == zero+layer)
00277 {
00278
00279 if(layer < max_radius-1)
00280 band_index.push_back(index);
00281
00282 int k = index / (width*height);
00283 int j = (index - k*width*height) / width;
00284 int i = index - (k*height+j)*width;
00285
00286
00287 for(unsigned int m=0; m<6; ++m)
00288 {
00289 int ii = i + neighbors[m][0];
00290 int jj = j + neighbors[m][1];
00291 int kk = k + neighbors[m][2];
00292
00293 if(ii >=0 && ii < (int)width && jj >= 0 && jj < (int)height && kk>=0 && kk<(int)depth)
00294 {
00295 int nindex = (kk*height+jj)*width+ii;
00296 if(phi[nindex] == zero)
00297 {
00298 phi[nindex] = zero+layer+1;
00299 inner_boundary[!current_list].push(nindex);
00300 }
00301 }
00302 }
00303 }
00304 }
00305
00306
00307 current_list = !current_list;
00308 ++layer;
00309 }
00310
00311
00312 current_list = 0;
00313 layer = 1;
00314 while(!outer_boundary[current_list].empty() && layer < max_radius)
00315 {
00316
00317 while(!outer_boundary[current_list].empty())
00318 {
00319
00320 int index = outer_boundary[current_list].front();
00321 outer_boundary[current_list].pop();
00322
00323
00324 if(phi[index] == zero-layer)
00325 {
00326
00327 if(layer < max_radius-1)
00328 band_index.push_back(index);
00329
00330 int k = index / (width*height);
00331 int j = (index - k*width*height) / width;
00332 int i = index - (k*height+j)*width;
00333
00334
00335 for(unsigned int m=0; m<6; ++m)
00336 {
00337 int ii = i + neighbors[m][0];
00338 int jj = j + neighbors[m][1];
00339 int kk = k + neighbors[m][2];
00340
00341 if(ii >=0 && ii < (int)width && jj >= 0 && jj < (int)height && kk>=0 && kk<(int)depth)
00342 {
00343 int nindex = (kk*height+jj)*width+ii;
00344 if(phi[nindex] == zero)
00345 {
00346 phi[nindex] = zero-layer-1;
00347 outer_boundary[!current_list].push(nindex);
00348 }
00349 }
00350 }
00351 }
00352 }
00353
00354
00355 current_list = !current_list;
00356 ++layer;
00357 }
00358 }
00359
00360
00361 template <typename DataType>
00362 void SmoothMask4D::march_cubes(const DataType* data,
00363 const unsigned char* mask,
00364 std::vector<int> &band_index,
00365 unsigned int width,
00366 unsigned int height,
00367 unsigned int depth)
00368 {
00369 try
00370 {
00371 float isovalue = 128;
00372
00373
00374 for(std::vector<int>::const_iterator it = band_index.begin(); it != band_index.end(); ++it)
00375 {
00376 int k = *it / (width*height);
00377 int j = (*it - k*width*height) / width;
00378 int i = *it - (k*height+j)*width;
00379
00380 if(i == (int)width-1 || j == (int)height-1 || k == (int)depth-1) continue;
00381
00382 GRIDCELL cell;
00383 cell.p[0].x = i+0; cell.p[0].y = j+0; cell.p[0].z = k+0; cell.val[0] = mask[((k+0)*height+(j+0))*width+(i+0)];
00384 cell.p[1].x = i+1; cell.p[1].y = j+0; cell.p[1].z = k+0; cell.val[1] = mask[((k+0)*height+(j+0))*width+(i+1)];
00385 cell.p[2].x = i+1; cell.p[2].y = j+1; cell.p[2].z = k+0; cell.val[2] = mask[((k+0)*height+(j+1))*width+(i+1)];
00386 cell.p[3].x = i+0; cell.p[3].y = j+1; cell.p[3].z = k+0; cell.val[3] = mask[((k+0)*height+(j+1))*width+(i+0)];
00387 cell.p[4].x = i+0; cell.p[4].y = j+0; cell.p[4].z = k+1; cell.val[4] = mask[((k+1)*height+(j+0))*width+(i+0)];
00388 cell.p[5].x = i+1; cell.p[5].y = j+0; cell.p[5].z = k+1; cell.val[5] = mask[((k+1)*height+(j+0))*width+(i+1)];
00389 cell.p[6].x = i+1; cell.p[6].y = j+1; cell.p[6].z = k+1; cell.val[6] = mask[((k+1)*height+(j+1))*width+(i+1)];
00390 cell.p[7].x = i+0; cell.p[7].y = j+1; cell.p[7].z = k+1; cell.val[7] = mask[((k+1)*height+(j+1))*width+(i+0)];
00391
00392 TRIANGLE triangles[5];
00393 int nbTriangles = 0;
00394
00395 nbTriangles = Polygonise(cell,isovalue,triangles);
00396
00397
00398 for(int m=0; m<nbTriangles; ++m)
00399 {
00400
00401 float x1 = triangles[m].p[0].x; float y1 = triangles[m].p[0].y; float z1 = triangles[m].p[0].z;
00402 float x2 = triangles[m].p[1].x; float y2 = triangles[m].p[1].y; float z2 = triangles[m].p[1].z;
00403 float x3 = triangles[m].p[2].x; float y3 = triangles[m].p[2].y; float z3 = triangles[m].p[2].z;
00404
00405
00406 float u[3] = {x2-x1, y2-y1, z2-z1};
00407 float v[3] = {x3-x2, y3-y2, z3-z2};
00408 float nx = + (u[1]*v[2] - v[1]*u[2]);
00409 float ny = - (u[0]*v[2] - v[0]*u[2]);
00410 float nz = + (u[0]*v[1] - v[0]*u[1]);
00411
00412
00413 if(!data)
00414 {
00415 glBegin(GL_TRIANGLES);
00416 glNormal3f(nx,ny,nz);
00417 glVertex3f(x1,y1,z1);
00418 glVertex3f(x2,y2,z2);
00419 glVertex3f(x3,y3,z3);
00420 glEnd();
00421 }
00422
00423
00424 else
00425 {
00426 glBegin(GL_TRIANGLES);
00427 for(unsigned int p=0; p<3; ++p)
00428 {
00429
00430 float x = triangles[m].p[p].x; float y = triangles[m].p[p].y; float z = triangles[m].p[p].z;
00431 int xp = (int)((x>0) ? x-1 : x);
00432 int xn = (int)((x<width-1) ? x+1 : x);
00433 int yp = (int)((y>0) ? y-1 : y);
00434 int yn = (int)((y<height-1) ? y+1 : y);
00435 int zp = (int)((z>0) ? z-1 : z);
00436 int zn = (int)((z<depth-1) ? z+1 : z);
00437
00438
00439 float gx = data[(((int)z)*height+((int)y))*width+(xn)] - data[(((int)z)*height+((int)y))*width+(xp)];
00440 float gy = data[(((int)z)*height+(yn))*width+((int)x)] - data[(((int)z)*height+(yp))*width+((int)x)];
00441 float gz = data[((zn)*height+((int)y))*width+((int)x)] - data[((zp)*height+((int)y))*width+((int)x)];
00442 float g = sqrt(gx*gx + gy*gy + gz*gz);
00443 gx = gx/g; gy = gy/g; gz = gz/g;
00444
00445
00446 float cos_theta = nx*gx + ny*gy + nz*gz;
00447 if(cos_theta < 0)
00448 {
00449 gx = -gx,
00450 gy = -gy,
00451 gz = -gz;
00452 }
00453
00454
00455 glNormal3f(gx,gy,gz);
00456 glVertex3f(x,y,z);
00457 }
00458 glEnd();
00459 }
00460 }
00461 }
00462 }
00463 catch(Error &e) { e.tag("IsoSurface::march_cubes"); throw e; }
00464 }
00465
00466 }
00467
00468 #endif