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

Isocontour.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 
00009 #include <stdio.h>
00010 #include <stdlib.h>
00011 #include <string.h>
00012 #include <math.h>
00013 #include "Inform.h"
00014 #include "utilities.h"
00015 
00016 #define ISOCONTOUR_INTERNAL 1
00017 #include "Isocontour.h"
00018 #include "VolumetricData.h"
00019 
00020 IsoContour::IsoContour(void) {}
00021 
00022 void IsoContour::clear(void) {
00023   numtriangles=0;
00024   v.clear();
00025   f.clear();
00026 }
00027 
00028 
00029 int IsoContour::compute(const VolumetricData *data, float isovalue, int step) {
00030   int x, y, z; 
00031   int tricount=0;
00032 
00033   vol=data;
00034 
00035   xax[0] = (float) (vol->xaxis[0] / ((float) (vol->xsize - 1)));
00036   xax[1] = (float) (vol->xaxis[1] / ((float) (vol->xsize - 1)));
00037   xax[2] = (float) (vol->xaxis[2] / ((float) (vol->xsize - 1)));
00038 
00039   yax[0] = (float) (vol->yaxis[0] / ((float) (vol->ysize - 1)));
00040   yax[1] = (float) (vol->yaxis[1] / ((float) (vol->ysize - 1)));
00041   yax[2] = (float) (vol->yaxis[2] / ((float) (vol->ysize - 1)));
00042 
00043   zax[0] = (float) (vol->zaxis[0] / ((float) (vol->zsize - 1)));
00044   zax[1] = (float) (vol->zaxis[1] / ((float) (vol->zsize - 1)));
00045   zax[2] = (float) (vol->zaxis[2] / ((float) (vol->zsize - 1)));
00046 
00047   for (z=0; z<(vol->zsize - step); z+=step) {
00048     for (y=0; y<(vol->ysize - step); y+=step) {
00049       for (x=0; x<(vol->xsize - step); x+=step) {
00050         tricount += DoCell(x, y, z, isovalue, step);
00051       }
00052     }
00053   }
00054 
00055   return 1;
00056 }
00057 
00058 
00059 
00060 int IsoContour::DoCell(int x, int y, int z, float isovalue, int step) {
00061   SQUARECELL gc;
00062   int addr, row, plane, rowstep, planestep, tricount;
00063   LINE tris[5];
00064 
00065   row = vol->xsize; 
00066   plane = vol->xsize * vol->ysize;
00067   addr = z*plane + y*row + x;
00068   rowstep = row*step;
00069   planestep = plane*step;
00070   
00071   gc.val[0] = vol->data[addr                             ];
00072   gc.val[1] = vol->data[addr + step                      ];
00073   gc.val[3] = vol->data[addr +        rowstep            ];
00074   gc.val[2] = vol->data[addr + step + rowstep            ];
00075   gc.val[4] = vol->data[addr +                  planestep];
00076   gc.val[5] = vol->data[addr + step +           planestep];
00077   gc.val[7] = vol->data[addr +        rowstep + planestep];
00078   gc.val[6] = vol->data[addr + step + rowstep + planestep];
00079 
00080   /*
00081      Determine the index into the edge table which
00082      tells us which vertices are inside of the surface
00083   */
00084   int cubeindex = 0;
00085   if (gc.val[0] < isovalue) cubeindex |= 1;
00086   if (gc.val[1] < isovalue) cubeindex |= 2;
00087   if (gc.val[2] < isovalue) cubeindex |= 4;
00088   if (gc.val[3] < isovalue) cubeindex |= 8;
00089   if (gc.val[4] < isovalue) cubeindex |= 16;
00090   if (gc.val[5] < isovalue) cubeindex |= 32;
00091   if (gc.val[6] < isovalue) cubeindex |= 64;
00092   if (gc.val[7] < isovalue) cubeindex |= 128;
00093 
00094   /* Cube is entirely in/out of the surface */
00095   if (edgeTable[cubeindex] == 0)
00096     return(0);
00097   gc.cubeindex = cubeindex;
00098 
00099   gc.p[0].x = (float) x;
00100   gc.p[0].y = (float) y;
00101 
00102   gc.p[1].x = (float) x + step;
00103   gc.p[1].y = (float) y;
00104 
00105   gc.p[2].x = (float) x + step;
00106   gc.p[2].y = (float) y + step;
00107   
00108   gc.p[3].x = (float) x;
00109   gc.p[3].y = (float) y + step;
00110 
00111   // calculate vertices and facets for this cube,
00112   // calculate normals by interpolating between the negated 
00113   // normalized volume gradients for the 8 reference voxels
00114   tricount = Polygonise(gc, isovalue, (LINE *) &tris);
00115 
00116   if (tricount > 0) {
00117     int i;
00118 
00119     for (i=0; i<tricount; i++) {
00120       float xx, yy, zz;
00121 
00122       int tritmp = numtriangles * 3;
00123       f.append3(tritmp, tritmp+1, tritmp+2);
00124       numtriangles++;
00125 
00126       // add new vertices and normals into vertex and normal lists
00127       xx = tris[i].p[0].x;
00128       yy = tris[i].p[0].y;
00129       v.append3((float) vol->origin[0] + xx * xax[0] + yy * yax[0] + zz * zax[0],
00130                 (float) vol->origin[1] + xx * xax[1] + yy * yax[1] + zz * zax[1],
00131                 (float) vol->origin[2] + xx * xax[2] + yy * yax[2] + zz * zax[2]);
00132       xx = tris[i].p[1].x;
00133       yy = tris[i].p[1].y;
00134       v.append3((float) vol->origin[0] + xx * xax[0] + yy * yax[0] + zz * zax[0],
00135                 (float) vol->origin[1] + xx * xax[1] + yy * yax[1] + zz * zax[1],
00136                 (float) vol->origin[2] + xx * xax[2] + yy * yax[2] + zz * zax[2]);
00137       xx = tris[i].p[2].x;
00138       yy = tris[i].p[2].y;
00139       v.append3((float) vol->origin[0] + xx * xax[0] + yy * yax[0] + zz * zax[0],
00140                 (float) vol->origin[1] + xx * xax[1] + yy * yax[1] + zz * zax[1],
00141                 (float) vol->origin[2] + xx * xax[2] + yy * yax[2] + zz * zax[2]);
00142     }
00143   }
00144 
00145   return tricount;
00146 }
00147 
00148 
00149 /*
00150    Given a grid cell and an isolevel, calculate the triangular
00151    facets required to represent the isocontour through the cell.
00152    Return the number of triangular facets, the array "triangles"
00153    will be loaded up with the vertices at most 5 triangular facets.
00154         0 will be returned if the grid cell is either totally above
00155    of totally below the isolevel.
00156    This code calculates vertex normals by interpolating the volume gradients.
00157 */
00158 int IsoContour::Polygonise(const SQUARECELL grid, const float isolevel, LINE *triangles) {
00159    int i,ntriang;
00160    int cubeindex = grid.cubeindex;
00161    XY vertlist[12];
00162 
00163    /* Find the vertices where the surface intersects the cube */
00164    if (edgeTable[cubeindex] & 1)
00165       VertexInterp(isolevel, grid, 0, 1, &vertlist[0]);
00166    if (edgeTable[cubeindex] & 2)
00167       VertexInterp(isolevel, grid, 1, 2, &vertlist[1]);
00168    if (edgeTable[cubeindex] & 4)
00169       VertexInterp(isolevel, grid, 2, 3, &vertlist[2]);
00170    if (edgeTable[cubeindex] & 8)
00171       VertexInterp(isolevel, grid, 3, 0, &vertlist[3]);
00172    if (edgeTable[cubeindex] & 16)
00173       VertexInterp(isolevel, grid, 4, 5, &vertlist[4]);
00174    if (edgeTable[cubeindex] & 32)
00175       VertexInterp(isolevel, grid, 5, 6, &vertlist[5]);
00176    if (edgeTable[cubeindex] & 64)
00177       VertexInterp(isolevel, grid, 6, 7, &vertlist[6]);
00178    if (edgeTable[cubeindex] & 128)
00179       VertexInterp(isolevel, grid, 7, 4, &vertlist[7]);
00180    if (edgeTable[cubeindex] & 256)
00181       VertexInterp(isolevel, grid, 0, 4, &vertlist[8]);
00182    if (edgeTable[cubeindex] & 512)
00183       VertexInterp(isolevel, grid, 1, 5, &vertlist[9]);
00184    if (edgeTable[cubeindex] & 1024)
00185       VertexInterp(isolevel, grid, 2, 6, &vertlist[10]);
00186    if (edgeTable[cubeindex] & 2048)
00187       VertexInterp(isolevel, grid, 3, 7, &vertlist[11]);
00188 
00189    /* Create the triangle */
00190    ntriang = 0;
00191    for (i=0; lineTable[cubeindex][i]!=-1;i+=3) {
00192      triangles[ntriang].p[0] = vertlist[lineTable[cubeindex][i  ]];
00193      triangles[ntriang].p[1] = vertlist[lineTable[cubeindex][i+1]];
00194      triangles[ntriang].p[2] = vertlist[lineTable[cubeindex][i+2]];
00195      ntriang++;
00196    }
00197 
00198    return ntriang;
00199 }
00200 
00201 
00202 /*
00203    Linearly interpolate the position where an isocontour cuts
00204    an edge between two vertices, each with their own scalar value,
00205    interpolating vertex position and vertex normal based on the 
00206    isovalue.
00207 */
00208 void IsoContour::VertexInterp(float isolevel, const SQUARECELL grid, int ind1, int ind2, XY * vert) {
00209   float mu;
00210 
00211   XY p1 = grid.p[ind1];
00212   XY p2 = grid.p[ind2];
00213   float valp1 = grid.val[ind1];
00214   float valp2 = grid.val[ind2];
00215 
00216   if (fabs(isolevel-valp1) < 0.00001) {
00217     *vert = grid.p[ind1];
00218     return;
00219   }
00220 
00221   if (fabs(isolevel-valp2) < 0.00001) {
00222     *vert = grid.p[ind2];
00223     return;
00224   }
00225 
00226   if (fabs(valp1-valp2) < 0.00001) {
00227     *vert = grid.p[ind1];
00228     return;
00229   }
00230 
00231   mu = (isolevel - valp1) / (valp2 - valp1);
00232 
00233   vert->x = p1.x + mu * (p2.x - p1.x);
00234   vert->y = p1.y + mu * (p2.y - p1.y);
00235 }
00236 

Generated on Fri Oct 4 02:44:09 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002