00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00027 #ifndef ISOSURFACE_H
00028 #define ISOSURFACE_H
00029
00030
00031 #ifdef __APPLE__
00032 #include <OpenGL/gl.h>
00033 #else
00034 #include <GL/gl.h>
00035 #endif
00036 #include <cmath>
00037 #include <error.h>
00038
00039
00040 namespace Display
00041 {
00076 class IsoSurface
00077 {
00078 public:
00079 IsoSurface();
00080 ~IsoSurface();
00081
00083 void display() const;
00084
00093 template <typename DataType>
00094 void update(const DataType* data,
00095 unsigned int width,
00096 unsigned int height,
00097 unsigned depth,
00098 float isovalue,
00099 bool from_bright_to_dark = true,
00100 bool use_gradient = true);
00101
00102 private:
00104 template <typename DataType>
00105 void march_cubes(const DataType* data,
00106 unsigned int width,
00107 unsigned int height,
00108 unsigned int depth,
00109 float isovalue,
00110 bool from_bright_to_dark,
00111 bool use_gradient);
00112
00115 struct XYZ { double x,y,z; };
00116 struct TRIANGLE { XYZ p[3]; };
00117 struct GRIDCELL { XYZ p[8]; double val[8]; };
00118 static XYZ VertexInterp(double isolevel, XYZ p1, XYZ p2, double valp1, double valp2);
00119 static int Polygonise(GRIDCELL grid,double isolevel,TRIANGLE *triangles);
00120
00121 private:
00122 GLuint m_meshid;
00123 };
00124
00125
00126 template <typename DataType>
00127 void IsoSurface::update(const DataType* data,
00128 unsigned int width,
00129 unsigned int height,
00130 unsigned depth,
00131 float isovalue,
00132 bool from_bright_to_dark,
00133 bool use_gradient)
00134 {
00135
00136 if(glIsList(m_meshid)) glDeleteLists(m_meshid,1);
00137 m_meshid = glGenLists(1);
00138 glNewList(m_meshid, GL_COMPILE);
00139
00140
00141 glMatrixMode(GL_MODELVIEW);
00142 glPushMatrix();
00143 glTranslatef(-0.5,-0.5,-0.5);
00144 glScalef(1.0f/width, 1.0f/height, 1.0f/depth);
00145
00146
00147 glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
00148 glEnable(GL_COLOR_MATERIAL);
00149
00150
00151 glEnable(GL_LIGHTING);
00152 glEnable(GL_LIGHT0);
00153
00154
00155 march_cubes(data, width, height, depth, isovalue, from_bright_to_dark, use_gradient);
00156
00157
00158 glDisable(GL_LIGHT0);
00159 glDisable(GL_LIGHTING);
00160 glDisable(GL_COLOR_MATERIAL);
00161
00162
00163 glPopMatrix();
00164 glEndList();
00165 }
00166
00167
00168
00169 template <typename DataType>
00170 void IsoSurface::march_cubes(const DataType* data,
00171 unsigned int width,
00172 unsigned int height,
00173 unsigned int depth,
00174 float isovalue,
00175 bool from_bright_to_dark,
00176 bool use_gradient)
00177 {
00178 try
00179 {
00180 for(unsigned int k=0; k<depth-1; ++k)
00181 for(unsigned int j=0; j<height-1; ++j)
00182 for(unsigned int i=0; i<width-1; ++i)
00183 {
00184 GRIDCELL cell;
00185 cell.p[0].x = i+0; cell.p[0].y = j+0; cell.p[0].z = k+0; cell.val[0] = data[((k+0)*height+(j+0))*width+(i+0)];
00186 cell.p[1].x = i+1; cell.p[1].y = j+0; cell.p[1].z = k+0; cell.val[1] = data[((k+0)*height+(j+0))*width+(i+1)];
00187 cell.p[2].x = i+1; cell.p[2].y = j+1; cell.p[2].z = k+0; cell.val[2] = data[((k+0)*height+(j+1))*width+(i+1)];
00188 cell.p[3].x = i+0; cell.p[3].y = j+1; cell.p[3].z = k+0; cell.val[3] = data[((k+0)*height+(j+1))*width+(i+0)];
00189 cell.p[4].x = i+0; cell.p[4].y = j+0; cell.p[4].z = k+1; cell.val[4] = data[((k+1)*height+(j+0))*width+(i+0)];
00190 cell.p[5].x = i+1; cell.p[5].y = j+0; cell.p[5].z = k+1; cell.val[5] = data[((k+1)*height+(j+0))*width+(i+1)];
00191 cell.p[6].x = i+1; cell.p[6].y = j+1; cell.p[6].z = k+1; cell.val[6] = data[((k+1)*height+(j+1))*width+(i+1)];
00192 cell.p[7].x = i+0; cell.p[7].y = j+1; cell.p[7].z = k+1; cell.val[7] = data[((k+1)*height+(j+1))*width+(i+0)];
00193
00194 TRIANGLE triangles[5];
00195 int nbTriangles = 0;
00196
00197 nbTriangles = Polygonise(cell,isovalue,triangles);
00198
00199
00200
00201 if(!use_gradient)
00202 {
00203 for(int m=0; m<nbTriangles; ++m)
00204 {
00205
00206 float x1 = triangles[m].p[0].x; float y1 = triangles[m].p[0].y; float z1 = triangles[m].p[0].z;
00207 float x2 = triangles[m].p[1].x; float y2 = triangles[m].p[1].y; float z2 = triangles[m].p[1].z;
00208 float x3 = triangles[m].p[2].x; float y3 = triangles[m].p[2].y; float z3 = triangles[m].p[2].z;
00209
00210
00211 float u[3] = {x2-x1, y2-y1, z2-z1};
00212 float v[3] = {x3-x2, y3-y2, z3-z2};
00213 float gx = + (u[1]*v[2] - v[1]*u[2]);
00214 float gy = - (u[0]*v[2] - v[0]*u[2]);
00215 float gz = + (u[0]*v[1] - v[0]*u[1]);
00216
00217
00218 glBegin(GL_TRIANGLES);
00219 glNormal3f(gx,gy,gz);
00220 glVertex3f(x1,y1,z1);
00221 glVertex3f(x2,y2,z2);
00222 glVertex3f(x3,y3,z3);
00223 glEnd();
00224 }
00225 }
00226 else
00227 {
00228 for(int m=0; m<nbTriangles; ++m)
00229 {
00230 glBegin(GL_TRIANGLES);
00231 for(unsigned int p=0; p<3; ++p)
00232 {
00233
00234 float x = triangles[m].p[p].x; float y = triangles[m].p[p].y; float z = triangles[m].p[p].z;
00235 int xp = (int)((x>0) ? x-1 : x);
00236 int xn = (int)((x<width-1) ? x+1 : x);
00237 int yp = (int)((y>0) ? y-1 : y);
00238 int yn = (int)((y<height-1) ? y+1 : y);
00239 int zp = (int)((z>0) ? z-1 : z);
00240 int zn = (int)((z<depth-1) ? z+1 : z);
00241
00242
00243 float gx = data[(((int)z)*height+((int)y))*width+(xn)] - data[(((int)z)*height+((int)y))*width+(xp)];
00244 float gy = data[(((int)z)*height+(yn))*width+((int)x)] - data[(((int)z)*height+(yp))*width+((int)x)];
00245 float gz = data[((zn)*height+((int)y))*width+((int)x)] - data[((zp)*height+((int)y))*width+((int)x)];
00246 float g = sqrt(gx*gx + gy*gy + gz*gz);
00247 gx = gx/g; gy = gy/g; gz = gz/g;
00248
00249
00250 if(from_bright_to_dark)
00251 gx = -gx,
00252 gy = -gy,
00253 gz = -gz;
00254
00255
00256 glNormal3f(gx,gy,gz);
00257 glVertex3f(x,y,z);
00258 }
00259 glEnd();
00260 }
00261 }
00262 }
00263 }
00264 catch(Error &e) { e.tag("IsoSurface::march_cubes"); throw e; }
00265 }
00266
00267 }
00268
00269 #endif