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