Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

Isosurface.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2019 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 #include <stdio.h>
00009 #include <stdlib.h>
00010 #include <string.h>
00011 #include <math.h>
00012 #include "Inform.h"
00013 #include "utilities.h"
00014 
00015 #define ISOSURFACE_INTERNAL 1
00016 #include "Isosurface.h"
00017 #include "VolumetricData.h"
00018 
00019 #define MIN(X,Y) (((X)<(Y))? (X) : (Y))
00020 #define MAX(X,Y) (((X)>(Y))? (X) : (Y))
00021 
00022 IsoSurface::IsoSurface(void) {}
00023 
00024 void IsoSurface::clear(void) {
00025   numtriangles=0;
00026   v.clear();
00027   n.clear();
00028   c.clear();
00029   f.clear();
00030 }
00031 
00032 
00033 int IsoSurface::compute(VolumetricData *data, float isovalue, int step) {
00034   int x, y, z; 
00035   int tricount=0;
00036 
00037   vol=data;
00038 
00039   // calculate cell axes
00040   vol->cell_axes(xax, yax, zax);
00041   vol->cell_dirs(xad, yad, zad);
00042 
00043   // flip normals if coordinate system is in the wrong handedness
00044   float vtmp[3];
00045   cross_prod(vtmp, xad, yad);
00046   if (dot_prod(vtmp, zad) < 0) {
00047     xad[0] *= -1;
00048     xad[1] *= -1;
00049     xad[2] *= -1;
00050     yad[0] *= -1;
00051     yad[1] *= -1;
00052     yad[2] *= -1;
00053     zad[0] *= -1;
00054     zad[1] *= -1;
00055     zad[2] *= -1;
00056   }
00057 
00058   // check that the grid is in the all-positive octant of the coordinate system
00059   int i;
00060   int axisposnorms = 1;
00061   for (i=0; i<3; i++) {
00062     if (xad[i] < 0 || yad[i] < 0 || zad[i] < 0)  
00063       axisposnorms = 0;
00064   }
00065 
00066   // check that the grid is axis-aligned
00067   if (xax[1] != 0.0f || xax[2] != 0.0f ||
00068       yax[0] != 0.0f || xax[2] != 0.0f ||
00069       yax[0] != 0.0f || xax[1] != 0.0f)
00070     axisposnorms = 0;
00071 
00072   if (axisposnorms) {
00073     tricount = DoGridPosNorms(isovalue, step);
00074   } else {
00075     // general case, any handedness, non-rectangular grids
00076     for (z=0; z<(vol->zsize - step); z+=step) {
00077       for (y=0; y<(vol->ysize - step); y+=step) {
00078         for (x=0; x<(vol->xsize - step); x+=step) {
00079           tricount += DoCellGeneral(x, y, z, isovalue, step);
00080         }
00081       }
00082     }
00083   }
00084 
00085   return 1;
00086 }
00087 
00088 
00089 //
00090 // Special case for axis-aligned grid with positive unit vector normals
00091 //
00092 int IsoSurface::DoGridPosNorms(float isovalue, int step) {
00093   GRIDCELL gc;
00094   TRIANGLE tris[5];
00095   long tricount=0;
00096   long globtricount=0;
00097 
00098   // precompute plane and row sizes to eliminate indirection
00099   long psz = vol->xsize * vol->ysize;
00100   long rsz = vol->xsize;
00101   long rowstep = rsz*step;
00102   long planestep = psz*step;
00103 
00104   // get direct access to gradient data for speed, and force 
00105   // generation of the gradients if they don't already exist
00106   const float *volgradient = vol->access_volume_gradient();
00107 
00108   int x, y, z;
00109   for (z=0; z<(vol->zsize - step); z+=step) {
00110     for (y=0; y<(vol->ysize - step); y+=step) {
00111       for (x=0; x<(vol->xsize - step); x+=step) {
00112         long addr = z*psz + y*rsz + x;
00113         gc.val[0] = vol->data[addr                             ];
00114         gc.val[1] = vol->data[addr + step                      ];
00115         gc.val[3] = vol->data[addr +        rowstep            ];
00116         gc.val[2] = vol->data[addr + step + rowstep            ];
00117         gc.val[4] = vol->data[addr +                  planestep];
00118         gc.val[5] = vol->data[addr + step +           planestep];
00119         gc.val[7] = vol->data[addr +        rowstep + planestep];
00120         gc.val[6] = vol->data[addr + step + rowstep + planestep];
00121 
00122         // Determine the index into the edge table which
00123         // tells us which vertices are inside of the surface
00124         int cubeindex = 0;
00125         if (gc.val[0] < isovalue) cubeindex |= 1;
00126         if (gc.val[1] < isovalue) cubeindex |= 2;
00127         if (gc.val[2] < isovalue) cubeindex |= 4;
00128         if (gc.val[3] < isovalue) cubeindex |= 8;
00129         if (gc.val[4] < isovalue) cubeindex |= 16;
00130         if (gc.val[5] < isovalue) cubeindex |= 32;
00131         if (gc.val[6] < isovalue) cubeindex |= 64;
00132         if (gc.val[7] < isovalue) cubeindex |= 128;
00133 
00134         // Cube is entirely in/out of the surface
00135         if (edgeTable[cubeindex] == 0)
00136           continue;
00137 
00138         gc.cubeindex = cubeindex;
00139 
00140         gc.p[0].x = (float) x;
00141         gc.p[0].y = (float) y;
00142         gc.p[0].z = (float) z;
00143         VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x, y, z, &gc.g[0].x)
00144 
00145         gc.p[1].x = (float) x + step;
00146         gc.p[1].y = (float) y;
00147         gc.p[1].z = (float) z;
00148         VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x + step, y, z, &gc.g[1].x)
00149 
00150         gc.p[3].x = (float) x;
00151         gc.p[3].y = (float) y + step;
00152         gc.p[3].z = (float) z;
00153         VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x, y + step, z, &gc.g[3].x)
00154 
00155         gc.p[2].x = (float) x + step;
00156         gc.p[2].y = (float) y + step;
00157         gc.p[2].z = (float) z;
00158         VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x + step, y + step, z, &gc.g[2].x)
00159 
00160         gc.p[4].x = (float) x;
00161         gc.p[4].y = (float) y;
00162         gc.p[4].z = (float) z + step;
00163         VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x, y, z + step, &gc.g[4].x)
00164 
00165         gc.p[5].x = (float) x + step;
00166         gc.p[5].y = (float) y;
00167         gc.p[5].z = (float) z + step;
00168         VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x + step, y, z + step, &gc.g[5].x)
00169 
00170         gc.p[7].x = (float) x;
00171         gc.p[7].y = (float) y + step;
00172         gc.p[7].z = (float) z + step;
00173         VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x, y + step, z + step, &gc.g[7].x)
00174 
00175         gc.p[6].x = (float) x + step;
00176         gc.p[6].y = (float) y + step;
00177         gc.p[6].z = (float) z + step;
00178         VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x + step, y + step, z + step, &gc.g[6].x)
00179 
00180         // calculate vertices and facets for this cube,
00181         // calculate normals by interpolating between the negated 
00182         // normalized volume gradients for the 8 reference voxels
00183         tricount = Polygonise(gc, isovalue, (TRIANGLE *) &tris);
00184         globtricount += tricount;
00185 
00186         // XXX this loop should avoid per-triangle append() operations
00187         //     in favor of directly writing to the output buffer, by 
00188         //     first calling the extend() method optimistically, and
00189         //     then truncating the output size to the actual triangle
00190         //     count that was emitted.  That would get the ResizeArray
00191         //     append/check methods entirely out of the innermost loops.
00192         int i;
00193         for (i=0; i<tricount; i++) {
00194           // XXX we need to truncate buffers larger than 32-bit ints since
00195           //     OpenGL etc won't like such large indices
00196           int tritmp = numtriangles * 3;
00197           f.append3(tritmp, tritmp+1, tritmp+2);
00198           numtriangles++;
00199 
00200           // add new vertices and normals into vertex and normal lists
00201           v.append3((float) vol->origin[0] + tris[i].p[0].x * xax[0],
00202                     (float) vol->origin[1] + tris[i].p[0].y * yax[1],
00203                     (float) vol->origin[2] + tris[i].p[0].z * zax[2]);
00204           n.append3(tris[i].n[0].x, tris[i].n[0].y, tris[i].n[0].z);
00205 
00206           v.append3((float) vol->origin[0] + tris[i].p[1].x * xax[0],
00207                     (float) vol->origin[1] + tris[i].p[1].y * yax[1],
00208                     (float) vol->origin[2] + tris[i].p[1].z * zax[2]);
00209           n.append3(tris[i].n[1].x, tris[i].n[1].y, tris[i].n[1].z);
00210 
00211           v.append3((float) vol->origin[0] + tris[i].p[2].x * xax[0],
00212                     (float) vol->origin[1] + tris[i].p[2].y * yax[1],
00213                     (float) vol->origin[2] + tris[i].p[2].z * zax[2]);
00214           n.append3(tris[i].n[2].x, tris[i].n[2].y, tris[i].n[2].z);
00215         }
00216       }
00217     }
00218   }
00219 
00220   return globtricount;
00221 }
00222 
00223 
00224 int IsoSurface::DoCellGeneral(int x, int y, int z, float isovalue, int step) {
00225   GRIDCELL gc;
00226   int tricount;
00227   TRIANGLE tris[5];
00228 
00229   // get direct access to gradient data for speed, and force 
00230   // generation of the gradients if they don't already exist
00231   const float *volgradient = vol->access_volume_gradient();
00232 
00233   // precompute plane/row sizes to eliminate indirection, promote int types
00234   long psz = vol->xsize * vol->ysize;
00235   long rsz = vol->xsize;
00236   long addr = z*psz + y*rsz + x;
00237   long rowstep = rsz*step;
00238   long planestep = psz*step;
00239   
00240   gc.val[0] = vol->data[addr                             ];
00241   gc.val[1] = vol->data[addr + step                      ];
00242   gc.val[3] = vol->data[addr +        rowstep            ];
00243   gc.val[2] = vol->data[addr + step + rowstep            ];
00244   gc.val[4] = vol->data[addr +                  planestep];
00245   gc.val[5] = vol->data[addr + step +           planestep];
00246   gc.val[7] = vol->data[addr +        rowstep + planestep];
00247   gc.val[6] = vol->data[addr + step + rowstep + planestep];
00248 
00249   // Determine the index into the edge table which
00250   // tells us which vertices are inside of the surface
00251   int cubeindex = 0;
00252   if (gc.val[0] < isovalue) cubeindex |= 1;
00253   if (gc.val[1] < isovalue) cubeindex |= 2;
00254   if (gc.val[2] < isovalue) cubeindex |= 4;
00255   if (gc.val[3] < isovalue) cubeindex |= 8;
00256   if (gc.val[4] < isovalue) cubeindex |= 16;
00257   if (gc.val[5] < isovalue) cubeindex |= 32;
00258   if (gc.val[6] < isovalue) cubeindex |= 64;
00259   if (gc.val[7] < isovalue) cubeindex |= 128;
00260 
00261   // Cube is entirely in/out of the surface
00262   if (edgeTable[cubeindex] == 0)
00263     return 0;
00264 
00265   gc.cubeindex = cubeindex;
00266 
00267   gc.p[0].x = (float) x;
00268   gc.p[0].y = (float) y;
00269   gc.p[0].z = (float) z;
00270   VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x, y, z, &gc.g[0].x)
00271 
00272   gc.p[1].x = (float) x + step;
00273   gc.p[1].y = (float) y;
00274   gc.p[1].z = (float) z;
00275   VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x + step, y, z, &gc.g[1].x)
00276 
00277   gc.p[3].x = (float) x;
00278   gc.p[3].y = (float) y + step;
00279   gc.p[3].z = (float) z;
00280   VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x, y + step, z, &gc.g[3].x)
00281 
00282   gc.p[2].x = (float) x + step;
00283   gc.p[2].y = (float) y + step;
00284   gc.p[2].z = (float) z;
00285   VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x + step, y + step, z, &gc.g[2].x)
00286 
00287   gc.p[4].x = (float) x;
00288   gc.p[4].y = (float) y;
00289   gc.p[4].z = (float) z + step;
00290   VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x, y, z + step, &gc.g[4].x)
00291 
00292   gc.p[5].x = (float) x + step;
00293   gc.p[5].y = (float) y;
00294   gc.p[5].z = (float) z + step;
00295   VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x + step, y, z + step, &gc.g[5].x)
00296 
00297   gc.p[7].x = (float) x;
00298   gc.p[7].y = (float) y + step;
00299   gc.p[7].z = (float) z + step;
00300   VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x, y + step, z + step, &gc.g[7].x)
00301 
00302   gc.p[6].x = (float) x + step;
00303   gc.p[6].y = (float) y + step;
00304   gc.p[6].z = (float) z + step;
00305   VOXEL_GRADIENT_FAST(volgradient, psz, rsz, x + step, y + step, z + step, &gc.g[6].x)
00306 
00307   // calculate vertices and facets for this cube,
00308   // calculate normals by interpolating between the negated 
00309   // normalized volume gradients for the 8 reference voxels
00310   tricount = Polygonise(gc, isovalue, (TRIANGLE *) &tris);
00311 
00312   // XXX this loop should avoid per-triangle append() operations
00313   //     in favor of directly writing to the output buffer, by 
00314   //     first calling the extend() method optimistically, and
00315   //     then truncating the output size to the actual triangle
00316   //     count that was emitted.  That would get the ResizeArray
00317   //     append/check methods entirely out of the innermost loops.
00318   int i;
00319   for (i=0; i<tricount; i++) {
00320     float xx, yy, zz;
00321     float xn, yn, zn;
00322 
00323     // XXX we need to truncate buffers larger than 32-bit ints since
00324     //     OpenGL etc won't like such large indices
00325     int tritmp = numtriangles * 3;
00326     f.append3(tritmp, tritmp+1, tritmp+2);
00327     numtriangles++;
00328 
00329     // add new vertices and normals into vertex and normal lists
00330     xx = tris[i].p[0].x;
00331     yy = tris[i].p[0].y;
00332     zz = tris[i].p[0].z;
00333 
00334     v.append3((float) vol->origin[0] + xx * xax[0] + yy * yax[0] + zz * zax[0],
00335               (float) vol->origin[1] + xx * xax[1] + yy * yax[1] + zz * zax[1],
00336               (float) vol->origin[2] + xx * xax[2] + yy * yax[2] + zz * zax[2]);
00337 
00338     xn = tris[i].n[0].x;
00339     yn = tris[i].n[0].y;
00340     zn = tris[i].n[0].z;
00341     n.append3((float) xn * xad[0] + yn * yad[0] + zn * zad[0],
00342               (float) xn * xad[1] + yn * yad[1] + zn * zad[1],
00343               (float) xn * xad[2] + yn * yad[2] + zn * zad[2]);
00344 
00345     xx = tris[i].p[1].x;
00346     yy = tris[i].p[1].y;
00347     zz = tris[i].p[1].z;
00348     v.append3((float) vol->origin[0] + xx * xax[0] + yy * yax[0] + zz * zax[0],
00349               (float) vol->origin[1] + xx * xax[1] + yy * yax[1] + zz * zax[1],
00350               (float) vol->origin[2] + xx * xax[2] + yy * yax[2] + zz * zax[2]);
00351 
00352     xn = tris[i].n[1].x;
00353     yn = tris[i].n[1].y;
00354     zn = tris[i].n[1].z;
00355     n.append3((float) xn * xad[0] + yn * yad[0] + zn * zad[0],
00356               (float) xn * xad[1] + yn * yad[1] + zn * zad[1],
00357               (float) xn * xad[2] + yn * yad[2] + zn * zad[2]);
00358 
00359     xx = tris[i].p[2].x;
00360     yy = tris[i].p[2].y;
00361     zz = tris[i].p[2].z;
00362     v.append3((float) vol->origin[0] + xx * xax[0] + yy * yax[0] + zz * zax[0],
00363               (float) vol->origin[1] + xx * xax[1] + yy * yax[1] + zz * zax[1],
00364               (float) vol->origin[2] + xx * xax[2] + yy * yax[2] + zz * zax[2]);
00365 
00366     xn = tris[i].n[2].x;
00367     yn = tris[i].n[2].y;
00368     zn = tris[i].n[2].z;
00369     n.append3((float) xn * xad[0] + yn * yad[0] + zn * zad[0],
00370               (float) xn * xad[1] + yn * yad[1] + zn * zad[1],
00371               (float) xn * xad[2] + yn * yad[2] + zn * zad[2]);
00372   }
00373 
00374   return tricount;
00375 }
00376 
00377 
00378 
00379 // normalize surface normals resulting from interpolation between 
00380 // unnormalized volume gradients
00381 void IsoSurface::normalize() {
00382   int i;
00383   int cnt = n.num();
00384   for (i=0; i<cnt; i+=3) {
00385     vec_normalize(&n[i]); 
00386   }  
00387 }
00388 
00389 
00390 // merge duplicated vertices detected by a simple windowed search
00391 int IsoSurface::vertexfusion(int offset, int len) {
00392   int i, j, newverts, oldverts, faceverts, matchcount;
00393 
00394   faceverts = f.num();
00395   oldverts = v.num() / 3; 
00396 
00397   // abort if we get an empty list
00398   if (!faceverts || !oldverts)
00399     return 0;
00400 
00401   int * vmap = new int[oldverts];
00402 
00403   vmap[0] = 0;
00404   newverts = 1;
00405   matchcount = 0;
00406 
00407   for (i=1; i<oldverts; i++) {
00408     int matchindex = -1;
00409     int start = ((newverts - offset) < 0)  ? 0        : (newverts - offset);
00410     int end   = ((start + len) > newverts) ? newverts : (start + len);
00411     int matched = 0;
00412     int vi = i * 3;
00413     for (j=start; j<end; j++) {
00414       int vj = j * 3;
00415       if (v[vi  ] == v[vj  ] && 
00416           v[vi+1] == v[vj+1] &&
00417           v[vi+2] == v[vj+2]) {
00418         matched = 1;
00419         matchindex = j;
00420         matchcount++;
00421         break;
00422       } 
00423     }
00424 
00425     if (matched) {
00426       vmap[i] = matchindex;
00427     } else {
00428       int vn = newverts * 3;
00429       v[vn    ] = v[vi    ];
00430       v[vn + 1] = v[vi + 1];
00431       v[vn + 2] = v[vi + 2];
00432       n[vn    ] = n[vi    ];
00433       n[vn + 1] = n[vi + 1];
00434       n[vn + 2] = n[vi + 2];
00435       vmap[i] = newverts;
00436       newverts++;
00437     } 
00438   }
00439 
00440 //  printf("Info) Vertex fusion: found %d shared vertices of %d, %d unique\n",
00441 //         matchcount, oldverts, newverts);
00442 
00443   // zap the old face, vertex, and normal arrays and replace with the new ones
00444   for (i=0; i<faceverts; i++) {
00445     f[i] = vmap[f[i]];
00446   }
00447   delete [] vmap;
00448 
00449   v.truncatelastn((oldverts - newverts) * 3);
00450   n.truncatelastn((oldverts - newverts) * 3);
00451 
00452   return 0;
00453 }
00454 
00455 
00456 /*
00457    Given a grid cell and an isolevel, calculate the triangular
00458    facets required to represent the isosurface through the cell.
00459    Return the number of triangular facets, the array "triangles"
00460    will be loaded up with the vertices at most 5 triangular facets.
00461         0 will be returned if the grid cell is either totally above
00462    of totally below the isolevel.
00463    This code calculates vertex normals by interpolating the volume gradients.
00464 */
00465 int IsoSurface::Polygonise(const GRIDCELL grid, const float isolevel, TRIANGLE *triangles) {
00466    int i,ntriang;
00467    int cubeindex = grid.cubeindex;
00468    XYZ vertlist[12];
00469    XYZ normlist[12];
00470 
00471    /* Find the vertices where the surface intersects the cube */
00472    if (edgeTable[cubeindex] & 1)
00473       VertexInterp(isolevel, grid, 0, 1, &vertlist[0], &normlist[0]);
00474    if (edgeTable[cubeindex] & 2)
00475       VertexInterp(isolevel, grid, 1, 2, &vertlist[1], &normlist[1]);
00476    if (edgeTable[cubeindex] & 4)
00477       VertexInterp(isolevel, grid, 2, 3, &vertlist[2], &normlist[2]);
00478    if (edgeTable[cubeindex] & 8)
00479       VertexInterp(isolevel, grid, 3, 0, &vertlist[3], &normlist[3]);
00480    if (edgeTable[cubeindex] & 16)
00481       VertexInterp(isolevel, grid, 4, 5, &vertlist[4], &normlist[4]);
00482    if (edgeTable[cubeindex] & 32)
00483       VertexInterp(isolevel, grid, 5, 6, &vertlist[5], &normlist[5]);
00484    if (edgeTable[cubeindex] & 64)
00485       VertexInterp(isolevel, grid, 6, 7, &vertlist[6], &normlist[6]);
00486    if (edgeTable[cubeindex] & 128)
00487       VertexInterp(isolevel, grid, 7, 4, &vertlist[7], &normlist[7]);
00488    if (edgeTable[cubeindex] & 256)
00489       VertexInterp(isolevel, grid, 0, 4, &vertlist[8], &normlist[8]);
00490    if (edgeTable[cubeindex] & 512)
00491       VertexInterp(isolevel, grid, 1, 5, &vertlist[9], &normlist[9]);
00492    if (edgeTable[cubeindex] & 1024)
00493       VertexInterp(isolevel, grid, 2, 6, &vertlist[10], &normlist[10]);
00494    if (edgeTable[cubeindex] & 2048)
00495       VertexInterp(isolevel, grid, 3, 7, &vertlist[11], &normlist[11]);
00496 
00497    /* Create the triangle */
00498    ntriang = 0;
00499    for (i=0;triTable[cubeindex][i]!=-1;i+=3) {
00500      triangles[ntriang].p[0] = vertlist[triTable[cubeindex][i  ]];
00501      triangles[ntriang].p[1] = vertlist[triTable[cubeindex][i+1]];
00502      triangles[ntriang].p[2] = vertlist[triTable[cubeindex][i+2]];
00503      triangles[ntriang].n[0] = normlist[triTable[cubeindex][i  ]];
00504      triangles[ntriang].n[1] = normlist[triTable[cubeindex][i+1]];
00505      triangles[ntriang].n[2] = normlist[triTable[cubeindex][i+2]];
00506      ntriang++;
00507    }
00508 
00509    return ntriang;
00510 }
00511 
00512 
00513 /*
00514    Linearly interpolate the position where an isosurface cuts
00515    an edge between two vertices, each with their own scalar value,
00516    interpolating vertex position and vertex normal based on the 
00517    isovalue.
00518 */
00519 void IsoSurface::VertexInterp(float isolevel, const GRIDCELL grid, int ind1, int ind2, XYZ * vert, XYZ * norm) {
00520   XYZ p1 = grid.p[ind1];
00521   XYZ p2 = grid.p[ind2];
00522   XYZ n1 = grid.g[ind1];
00523   XYZ n2 = grid.g[ind2];
00524   float valp1 = grid.val[ind1];
00525   float valp2 = grid.val[ind2];
00526   float isodiffp1   = isolevel - valp1;
00527   float diffvalp2p1 = valp2 - valp1;
00528   float mu = 0.0f;
00529 
00530   // if the difference between vertex values is zero or nearly
00531   // zero, we can get an IEEE NAN for mu.  We must either avoid this
00532   // by testing the denominator beforehand, by coping with the resulting
00533   // NAN value after the fact.  The only important thing is that mu be
00534   // assigned a value between zero and one.
00535   
00536 #if 0
00537   if (fabsf(isodiffp1) < 0.00001) {
00538     *vert = p1;
00539     *norm = n1;
00540     return;
00541   }
00542 
00543   if (fabsf(isolevel-valp2) < 0.00001) {
00544     *vert = p2;
00545     *norm = n2;
00546     return;
00547   }
00548 
00549   if (fabsf(diffvalp2p1) < 0.00001) {
00550     *vert = p1;
00551     *norm = n1;
00552     return;
00553   }
00554 #endif
00555 
00556   if (fabsf(diffvalp2p1) > 0.0f) 
00557     mu = isodiffp1 / diffvalp2p1;
00558 
00559 #if 0
00560   if (mu > 1.0f)
00561     mu=1.0f;
00562 
00563   if (mu < 0.0f)
00564     mu=0.0f;
00565 #endif
00566 
00567   vert->x = p1.x + mu * (p2.x - p1.x);
00568   vert->y = p1.y + mu * (p2.y - p1.y);
00569   vert->z = p1.z + mu * (p2.z - p1.z);
00570 
00571   norm->x = n1.x + mu * (n2.x - n1.x);
00572   norm->y = n1.y + mu * (n2.y - n1.y);
00573   norm->z = n1.z + mu * (n2.z - n1.z);
00574 }
00575 
00576 
00582 int IsoSurface::set_color_rgb3fv(const float *rgb) {
00583   int i;
00584   int numcols = v.num();
00585 
00586   // Pre-extend our output buffers so we don't have to call append().
00587   c.extend(numcols);
00588   for (i=0; i<numcols; i+=3) {
00589     c[i  ] = rgb[0];
00590     c[i+1] = rgb[1];
00591     c[i+2] = rgb[2];
00592   }
00593   c.set_size(i);
00594 
00595   return 0;
00596 }
00597 
00598 
00602 int IsoSurface::set_color_voltex_rgb3fv(const float *voltex) {
00603   int i;
00604   int numverts = v.num() / 3;
00605 
00606   // precompute plane and row sizes to eliminate indirection
00607   long psz = vol->xsize * vol->ysize;
00608   long rsz = vol->xsize; 
00609  
00610   float xinv = 1.0f / xax[0];
00611   float yinv = 1.0f / yax[1];
00612   float zinv = 1.0f / zax[2];
00613   int xs = vol->xsize - 1;
00614   int ys = vol->ysize - 1;
00615   int zs = vol->zsize - 1;
00616 
00617   for (i=0; i<numverts; i++) {
00618     long ind = i*3;
00619     float vx = float((v[ind    ] - vol->origin[0]) * xinv);
00620     float vy = float((v[ind + 1] - vol->origin[1]) * yinv);
00621     float vz = float((v[ind + 2] - vol->origin[2]) * zinv);
00622 
00623     int x = MIN(MAX(((int) vx), 0), xs);
00624     int y = MIN(MAX(((int) vy), 0), ys);
00625     int z = MIN(MAX(((int) vz), 0), zs);
00626 
00627     long caddr = (z*psz + y*rsz + x)*3L;
00628 
00629 #if 0
00630     // non-interpolated texture lookup
00631     float rgb[3];
00632     vec_copy(rgb, &voltex[caddr]);
00633 
00634     // normalize colors to eliminate gridding artifacts for the time being
00635     vec_normalize(rgb);
00636 #else
00637     float mux = (vx - x);
00638     float muy = (vy - y);
00639     float muz = (vz - z);
00640 
00641     float c0[3], c1[3], c2[3], c3[3]; 
00642 
00643     long caddrp1 = (x < xs) ? caddr+3 : caddr;
00644     c0[0] = (1-mux)*voltex[caddr    ] + mux*voltex[caddrp1    ];
00645     c0[1] = (1-mux)*voltex[caddr + 1] + mux*voltex[caddrp1 + 1];
00646     c0[2] = (1-mux)*voltex[caddr + 2] + mux*voltex[caddrp1 + 2];
00647 
00648     long yinc = (y < ys) ? rsz*3L : 0; 
00649     caddr += yinc;
00650     caddrp1 += yinc;
00651     c1[0] = (1-mux)*voltex[caddr    ] + mux*voltex[caddrp1    ];
00652     c1[1] = (1-mux)*voltex[caddr + 1] + mux*voltex[caddrp1 + 1];
00653     c1[2] = (1-mux)*voltex[caddr + 2] + mux*voltex[caddrp1 + 2];
00654     
00655     long zinc = (z < zs) ? psz*3L : 0; 
00656     caddr += zinc;
00657     caddrp1 += zinc;
00658     c3[0] = (1-mux)*voltex[caddr    ] + mux*voltex[caddrp1    ];
00659     c3[1] = (1-mux)*voltex[caddr + 1] + mux*voltex[caddrp1 + 1];
00660     c3[2] = (1-mux)*voltex[caddr + 2] + mux*voltex[caddrp1 + 2];
00661 
00662     caddr -= yinc;
00663     caddrp1 -= yinc;
00664     c2[0] = (1-mux)*voltex[caddr    ] + mux*voltex[caddrp1    ];
00665     c2[1] = (1-mux)*voltex[caddr + 1] + mux*voltex[caddrp1 + 1];
00666     c2[2] = (1-mux)*voltex[caddr + 2] + mux*voltex[caddrp1 + 2];
00667 
00668     float cz0[3];
00669     cz0[0] = (1-muy)*c0[0] + muy*c1[0];
00670     cz0[1] = (1-muy)*c0[1] + muy*c1[1];
00671     cz0[2] = (1-muy)*c0[2] + muy*c1[2];
00672 
00673     float cz1[3];
00674     cz1[0] = (1-muy)*c2[0] + muy*c3[0];
00675     cz1[1] = (1-muy)*c2[1] + muy*c3[1];
00676     cz1[2] = (1-muy)*c2[2] + muy*c3[2];
00677    
00678     float rgb[3];
00679     rgb[0] = (1-muz)*cz0[0] + muz*cz1[0];
00680     rgb[1] = (1-muz)*cz0[1] + muz*cz1[1];
00681     rgb[2] = (1-muz)*cz0[2] + muz*cz1[2];
00682 #endif
00683 
00684     c.append3(&rgb[0]);
00685   }
00686 
00687   return 0;
00688 }
00689 
00690 
00691 
00692 

Generated on Thu Mar 28 02:43:20 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002