/* vtmcP.c */ /* * Vis5D system for visualizing five dimensional gridded data sets. * Copyright (C) 1990 - 2000 Bill Hibbard, Johan Kellum, Brian Paul, * Dave Santek, and Andre Battaiola. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * As a special exception to the terms of the GNU General Public * License, you are permitted to link Vis5D with (and distribute the * resulting source and executables) the LUI library (copyright by * Stellar Computer Inc. and licensed for distribution with Vis5D), * the McIDAS library, and/or the NetCDF library, where those * libraries are governed by the terms of their own licenses. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA * */ #include "../config.h" /* ---------|---------|---------|---------|---------|---------|---------|-| PROGRAM TO CALCULATE MARCHCUBE CUBES NEW COMMENTS AND CHANGES: Andre Luiz Battaiola SSEC/INPE MODIFIED BY Brian Paul TO BE RE-ENTRANT AND ALLOCATE ALL MEMORY DYNAMICALLY. BUG FIXES: Checks for invalid polygon when arrays filled. BEP 2-13-92 -----------------------------------------------------------------------| */ #include #include #include #ifdef LEVELTYPES /* This is a kludge! @@ It was NOT possible to include globals.h here! */ /* API - note */ extern float Height[]; extern int LevelType; # define PRESSURE_LEVELS 100 # define HEIGHT_LEVELS 105 #endif #include "etableP.h" #include "memory.h" #ifndef TRUE #define TRUE ( 1 == 1 ) #define FALSE ( 1 == 0 ) #endif #ifndef min #define min(a,b) ((a < b) ? a : b) #endif #ifndef max #define max(a,b) ((a > b) ? a : b) #endif #define mod(a,b) (a==b ? 0 : a) #define rabs(a) ((a >= 0.) ? a : -a) #define BIG_NEG -2e+9 #define EPS_0 1.0e-5 #define EPS_1 1.0 - EPS_0 #define INVALID_VALUE 1.0e30 #define exist_polygon_in_cube(nc) \ ((ptFLAG[nc] != 0) & (ptFLAG[nc] != 0xFF)) /* API - global to all Vis5D contexts, and DEBUG implies a single work task */ #ifdef DEBUG static FILE *output; /* Pointer to the output file */ static void statistics(); #endif /* ----- Routine to allocate memory space with check */ static void *xalloc( Context ctx, int size, int type ) { return (void *) allocate_type( ctx, size, type); } /* ----- Routine to deallocate memory space with check */ static void xfree( Context ctx, void *pt ) { if (pt != NULL) deallocate(ctx, pt, -1); } /***********************************************************/ /**** FLAGS ****/ /***********************************************************/ /* ptGRID(x,y,z) -> *( ptGRID + y + x*ydim + z*xdim_x_ydim ) */ #define value_plane_z0_00 ptAUX[ pcube[ii] ] #define value_plane_z0_10 ptAUX[ pcube[ii] + ydim ] #define value_plane_z0_01 ptAUX[ pcube[ii] + 1 ] #define value_plane_z0_11 ptAUX[ pcube[ii] + ydim + 1 ] #define value_plane_z1_00 ptAUX[ pcube[ii] + xdim_x_ydim ] #define value_plane_z1_10 ptAUX[ pcube[ii] + ydim + xdim_x_ydim ] #define value_plane_z1_01 ptAUX[ pcube[ii] + 1 + xdim_x_ydim ] #define value_plane_z1_11 ptAUX[ pcube[ii] + 1 + ydim + xdim_x_ydim ] #define exist_cube (cb > -1 && cb < num_cubes) #define cZp cb #define cZn cb #define cYp cb #define cYn cb #define cXp cb #define cXn cb #define get_xyz_cube() \ { iz = ii / num_cubes_xy; \ ix = (int)((ii - (iz * num_cubes_xy)) / num_cubes_y); \ iy = ii - (iz * num_cubes_xy) - (ix * num_cubes_y); \ } #define get_cZp() \ cb = (iz < (zdim - 1)) ? ii + num_cubes_xy : -1 #define get_cZn() \ cb = (iz > 0) ? ii - num_cubes_xy : -1 #define get_cXp() \ cb = (ix < (xdim - 1)) ? ii + num_cubes_y : -1 #define get_cXn() \ cb = (ix > 0) ? ii - num_cubes_y : -1 #define get_cYp() \ cb = (iy < (ydim - 1)) ? ii + 1 : -1 #define get_cYn() \ cb = (iy > 0) ? ii - 1 : -1 /* ---------|---------|---------|---------|---------|---------|---------|-| Routine to Calculate FLAGS used in Marching Cubes Program Author: Andre Luiz Battaiola SSEC/INPE -----------------------------------------------------------------------| */ static int flags( float *ptGRID, int xdim, int ydim, int zdim, int ptFLAG[], int ptAUX[], int pcube[], double isovalue ) { int ii, jj, ix, iy, iz, cb, SF, bcase; int num_cubes, num_cubes_xy, num_cubes_y; int xdim_x_ydim = xdim*ydim; int xdim_x_ydim_x_zdim = xdim_x_ydim*zdim; int npolygons; num_cubes_y = ydim-1; num_cubes_xy = (xdim-1) * num_cubes_y; num_cubes = (zdim-1) * num_cubes_xy; /* ************* Big Attention ************* pcube must have the dimension "num_cubes+1" because in terms of eficiency the best loop to calculate "pcube" will use more one component to pcube. */ /* Calculate the Flag Value of each Cube */ /* In order to simplify the Flags calculations, "pcube" will be used to store the number of the first vertex of each cube */ ii = 0; pcube[0] = 0; for (iz=0; iz<(zdim-1); iz++) { for (ix=0; ix<(xdim-1); ix++) { cb = pcube[ii]; for (iy=1; iy<(ydim-1); iy++) /* Vectorized */ pcube[ii+iy] = cb+iy; ii += ydim-1; pcube[ii] = pcube[ii-1]+2; } pcube[ii] += ydim; } /* Vectorized */ for (ii = 0; ii < xdim_x_ydim_x_zdim; ii++) if (ptGRID[ii] >= INVALID_VALUE) ptAUX[ii] = 0x1001; else if (ptGRID[ii] >= isovalue) ptAUX[ii] = 1; else ptAUX[ii] = 0; /* Vectorized */ for (ii = 0; ii < num_cubes; ii++) ptFLAG[ii] = ((value_plane_z0_00 ) | (value_plane_z0_10 << 1) | (value_plane_z0_01 << 2) | (value_plane_z0_11 << 3) | (value_plane_z1_00 << 4) | (value_plane_z1_10 << 5) | (value_plane_z1_01 << 6) | (value_plane_z1_11 << 7)); /* After this Point it is not more used pcube */ /* Analyse Special Cases in FLAG */ ii = npolygons = 0; while ( TRUE ) { for (; ii < num_cubes; ii++) if (exist_polygon_in_cube(ii) && (ptFLAG[ii] < MAX_FLAG_NUM)) break; if ( ii == num_cubes ) break; bcase = pol_edges[ptFLAG[ii]][0]; if (bcase == 0xE6 || bcase == 0xF9) { get_xyz_cube(); /* == Z+ == */ if ((ptFLAG[ii] & 0xF0) == 0x90 || (ptFLAG[ii] & 0xF0) == 0x60) get_cZp(); /* == Z- == */ else if ((ptFLAG[ii] & 0x0F) == 0x09 || (ptFLAG[ii] & 0x0F) == 0x06) get_cZn(); /* == Y+ == */ else if ((ptFLAG[ii] & 0xCC) == 0x84 || (ptFLAG[ii] & 0xCC) == 0x48) get_cYp(); /* == Y- == */ else if ((ptFLAG[ii] & 0x33) == 0x21 || (ptFLAG[ii] & 0x33) == 0x12) get_cYn(); /* == X+ == */ else if ((ptFLAG[ii] & 0xAA) == 0x82 || (ptFLAG[ii] & 0xAA) == 0x28) get_cXp(); /* == X- == */ else if ((ptFLAG[ii] & 0x55) == 0x41 || (ptFLAG[ii] & 0x55) == 0x14) get_cXn(); /* == Map Special Case == */ if (exist_cube && ptFLAG[cb]<316) /*changed by BEP on 7-20-92*/ { bcase = pol_edges[ptFLAG[cb]][0]; if (bcase == 0x06 || bcase == 0x16 || bcase == 0x19 || bcase == 0x1E || bcase == 0x3C || bcase == 0x69) ptFLAG[ii] = sp_cases[ptFLAG[ii]]; } } else if (bcase == 0xE9) { get_xyz_cube(); if (ptFLAG[ii] == 0x6B) SF = SF_6B; else if (ptFLAG[ii] == 0x6D) SF = SF_6D; else if (ptFLAG[ii] == 0x79) SF = SF_79; else if (ptFLAG[ii] == 0x97) SF = SF_97; else if (ptFLAG[ii] == 0x9E) SF = SF_9E; else if (ptFLAG[ii] == 0xB6) SF = SF_B6; else if (ptFLAG[ii] == 0xD6) SF = SF_D6; else if (ptFLAG[ii] == 0xE9) SF = SF_E9; for (jj=0; jj<3; jj++) { if (case_E9[jj+SF] == Zp) get_cZp(); else if (case_E9[jj+SF] == Zn) get_cZn(); else if (case_E9[jj+SF] == Yp) get_cYp(); else if (case_E9[jj+SF] == Yn) get_cYn(); else if (case_E9[jj+SF] == Xp) get_cXp(); else if (case_E9[jj+SF] == Xn) get_cXn(); /* changed: if (exist_cube) to: */ if (exist_cube && ptFLAG[cb]<316) /*changed by BEP on 7-20-92*/ { bcase = pol_edges[ptFLAG[cb]][0]; if (bcase == 0x06 || bcase == 0x16 || bcase == 0x19 || bcase == 0x1E || bcase == 0x3C || bcase == 0x69) { ptFLAG[ii] = sp_cases[ptFLAG[ii]] + case_E9[jj+SF+3]; break; } } } } /* Calculate the Number of Generated Triangles and Polygons */ npolygons += pol_edges[ptFLAG[ii]][1]; ii++; } /* npolygons2 = 2*npolygons; */ return npolygons; } /********************************************************/ /*** MARCHING ****/ /********************************************************/ #define MASK 0x0F #define P_F_V(v,t) Pol_f_Vert[v*9 + t] #define V_f_P Vert_f_Pol #define num_polygons_in_cube(nc) pol_edges[ptFLAG[nc]][1] #ifdef DEBUG /* ----- Print Macro Routines - They are acessed only in case of the debug option has been setted */ #define print_new_cube(ncube,caseH) \ { if (caseH == 0x00) caseA = 0; \ else if (caseH == 0x01) caseA = 1; \ else if (caseH == 0x03) caseA = 2; \ else if (caseH == 0x06) caseA = 3; \ else if (caseH == 0x18) caseA = 4; \ else if (caseH == 0x07) caseA = 5; \ else if (caseH == 0x19) caseA = 6; \ else if (caseH == 0x16) caseA = 7; \ else if (caseH == 0x0f) caseA = 8; \ else if (caseH == 0x17) caseA = 9; \ else if (caseH == 0x3c) caseA = 10; \ else if (caseH == 0x1b) caseA = 11; \ else if (caseH == 0x1e) caseA = 12; \ else if (caseH == 0x69) caseA = 13; \ else if (caseH == 0x27) caseA = 14; \ fprintf (output,"\nCube Number: %d\n",ncube+1); \ fprintf (output, \ "Case: Program = %x, Article = %d\n", \ caseH,caseA); \ } #define print_no_polygon() \ fprintf (output,"-- no polygons generated --\n") #define print_new_polygon( m ) \ fprintf (output,"Polygon Number: %d\n",m) #define print_vertex( m ) \ fprintf(output,"%f %f %f\n",VX[m],VY[m],VZ[m]) #endif /* ---------|---------|---------|---------|---------|---------|---------|---------| */ /* v v ----- For one of the axis, this macro calculates the point that the curve will cross the line (one side of the cube) that join two points of this axis. Example: Suppose: a = x1, b = x2, dx = (x2 - x1) = 1, dn = (node_value_x2 - node_value_x1) so: xcross = ((isovalue - node_value_x1) * dx / dn) + x1 xcross = ((isovalue - node_value_x1) / dn) + x1 */ #define calcNode(i1,vi2,vi1) \ { nodeDiff = vi2 - vi1; \ cp = ( ( isovalue - vi1 ) / nodeDiff ) + i1; \ } /* ^ ^ ----- end calcNode */ /* printf("Old cp = %f ", (isovalue-node_below)/(node_above-node_below) + gridbox); \ printf("New cp = %f\n", cp); \ */ #ifdef LEVELTYPES #define LOGcalcNode(gridbox, node_above, node_below) \ { float Pbelow, Pabove, lnPwanted, Pwanted, f; \ \ if (node_above == node_below) \ f = 1.0; \ else \ f = (isovalue - node_below) / (node_above - node_below); \ \ Pbelow = Height[gridbox]; \ Pabove = Height[gridbox+1]; \ if (Pbelow == Pabove) \ cp = gridbox; \ else { \ lnPwanted = log(Pbelow) - f*(log(Pbelow) - log(Pabove)); \ Pwanted = exp(lnPwanted); \ cp = (Pbelow - Pwanted) / (Pbelow - Pabove) + gridbox; \ } \ } /* ^ ^ ----- end LOGcalcNode */ #endif #define valid_cube(nc) (ptFLAG[ncube] < MAX_FLAG_NUM) #define PlaneX(n,x,y) ixPlane[n] + x*ydim + y #define PlaneY(n,x,y) iyPlane[n] + y*xdim + x #define PlaneZ(n,y) izPlane[n] + y #define exist_edge_1 pol_edges[ptFLAG[ncube]][3] & 0x0002 #define exist_edge_2 pol_edges[ptFLAG[ncube]][3] & 0x0004 #define exist_edge_3 pol_edges[ptFLAG[ncube]][3] & 0x0008 #define exist_edge_4 pol_edges[ptFLAG[ncube]][3] & 0x0010 #define exist_edge_5 pol_edges[ptFLAG[ncube]][3] & 0x0020 #define exist_edge_6 pol_edges[ptFLAG[ncube]][3] & 0x0040 #define exist_edge_7 pol_edges[ptFLAG[ncube]][3] & 0x0080 #define exist_edge_8 pol_edges[ptFLAG[ncube]][3] & 0x0100 #define exist_edge_9 pol_edges[ptFLAG[ncube]][3] & 0x0200 #define exist_edge_A pol_edges[ptFLAG[ncube]][3] & 0x0400 #define exist_edge_B pol_edges[ptFLAG[ncube]][3] & 0x0800 #define exist_edge_C pol_edges[ptFLAG[ncube]][3] & 0x1000 #define vnode0 *(pt) #define vnode1 *(pt + ydim) #define vnode2 *(pt + 1) #define vnode3 *(pt + ydim + 1) #define vnode4 *(pt + xdim_x_ydim) #define vnode5 *(pt + ydim + xdim_x_ydim) #define vnode6 *(pt + 1 + xdim_x_ydim) #define vnode7 *(pt + 1 + ydim + xdim_x_ydim) /* ************* Big Attention ************* The "ifs" in the "find_vertex" macro was defined supposing that the cubes are processed first in "y" axes direction, after "x" axes direction, and finally in "z" direction. If that order is changed, it's necessary to redefine the "ifs". Remind that "Fortran" alignes arrays by columns and "C" alignes by rows. */ #define nvet nvertex #define find_vertex() \ { if (exist_edge_1) /* cube vertex 0-1 */ \ { if (iz || iy) calc_edge[1] = *(PlaneX(bellow,ix,iy)); \ else \ { calcNode(ix,vnode1,vnode0); calc_edge[1] = nvet; \ VX[nvet] = cp; VY[nvet] = iy; VZ[nvet] = iz; \ nvet++; \ } \ } \ if (exist_edge_2) /* cube vertex 0-2 */ \ { if (iz || ix) calc_edge[2] = *(PlaneY(bellow,ix,iy)); \ else \ { calcNode(iy,vnode2,vnode0); calc_edge[2] = nvet; \ VX[nvet] = ix; VY[nvet] = cp; VZ[nvet] = iz; \ nvet++; \ } \ } \ if (exist_edge_3) /* cube vertex 0-4 */ \ { if (ix || iy) calc_edge[3] = *(PlaneZ(rear,iy)); \ else \ { calcNode(iz,vnode4,vnode0); calc_edge[3] = nvet; \ VX[nvet] = ix; VY[nvet] = iy; VZ[nvet] = cp; \ nvet++; \ } \ } \ if (exist_edge_4) /* cube vertex 1-3 */ \ { if (iz) calc_edge[4] = *(PlaneY(bellow,ix+1,iy)); \ else \ { calcNode(iy,vnode3,vnode1); calc_edge[4] = nvet; \ VX[nvet] = ix+1; VY[nvet] = cp; VZ[nvet] = iz; \ *(PlaneY(bellow,ix+1,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_5) /* cube vertex 1-5 */ \ { if (iy) calc_edge[5] = *(PlaneZ(front,iy)); \ else \ { calcNode(iz,vnode5,vnode1); calc_edge[5] = nvet; \ VX[nvet] = ix+1; VY[nvet] = iy; VZ[nvet] = cp; \ *(PlaneZ(front,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_6) /* cube vertex 2-3 */ \ { if (iz) calc_edge[6] = *(PlaneX(bellow,ix,iy+1)); \ else \ { calcNode(ix,vnode3,vnode2); calc_edge[6] = nvet; \ VX[nvet] = cp; VY[nvet] = iy+1; VZ[nvet] = iz; \ *(PlaneX(bellow,ix,iy+1)) = nvet; nvet++; \ } \ } \ if (exist_edge_7) /* cube vertex 2-6 */ \ { if (ix) calc_edge[7] = *(PlaneZ(rear,iy+1)); \ else \ { calcNode(iz,vnode6,vnode2); calc_edge[7] = nvet; \ VX[nvet] = ix; VY[nvet] = iy+1; VZ[nvet] = cp; \ *(PlaneZ(rear,iy+1)) = nvet; nvet++; \ } \ } \ if (exist_edge_8) /* cube vertex 3-7 */ \ { calcNode(iz,vnode7,vnode3); calc_edge[8] = nvet; \ VX[nvet] = ix+1; VY[nvet] = iy+1; VZ[nvet] = cp; \ *(PlaneZ(front,iy+1)) = nvet; nvet++; \ } \ if (exist_edge_9) /* cube vertex 4-5 */ \ { if (iy) calc_edge[9] = *(PlaneX(above,ix,iy)); \ else \ { calcNode(ix,vnode5,vnode4); calc_edge[9] = nvet; \ VX[nvet] = cp; VY[nvet] = iy; VZ[nvet] = iz+1; \ *(PlaneX(above,ix,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_A) /* cube vertex 4-6 */ \ { if (ix) calc_edge[10] = *(PlaneY(above,ix,iy)); \ else \ { calcNode(iy,vnode6,vnode4); calc_edge[10] = nvet; \ VX[nvet] = ix; VY[nvet] = cp; VZ[nvet] = iz+1; \ *(PlaneY(above,ix,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_B) /* cube vertex 5-7 */ \ { calcNode(iy,vnode7,vnode5); calc_edge[11] = nvet; \ VX[nvet] = ix+1; VY[nvet] = cp; VZ[nvet] = iz+1; \ *(PlaneY(above,ix+1,iy)) = nvet; nvet++; \ } \ if (exist_edge_C) /* cube vertex 6-7 */ \ { calcNode(ix,vnode7,vnode6); calc_edge[12] = nvet; \ VX[nvet] = cp; VY[nvet] = iy+1; VZ[nvet] = iz+1; \ *(PlaneX(above,ix,iy+1)) = nvet; nvet++; \ } \ } #ifdef LEVELTYPES /* LOGfind_vertex is equal to find_vertex except for the macro's called. When interpolation in the Z-coordinate must be performed, LOGcalcNode is used. */ #define LOGfind_vertex() \ { if (exist_edge_1) /* cube vertex 0-1 */ \ { if (iz || iy) calc_edge[1] = *(PlaneX(bellow,ix,iy)); \ else \ { calcNode(ix,vnode1,vnode0); calc_edge[1] = nvet; \ VX[nvet] = cp; VY[nvet] = iy; VZ[nvet] = iz; \ nvet++; \ } \ } \ if (exist_edge_2) /* cube vertex 0-2 */ \ { if (iz || ix) calc_edge[2] = *(PlaneY(bellow,ix,iy)); \ else \ { calcNode(iy,vnode2,vnode0); calc_edge[2] = nvet; \ VX[nvet] = ix; VY[nvet] = cp; VZ[nvet] = iz; \ nvet++; \ } \ } \ if (exist_edge_3) /* cube vertex 0-4 */ \ { if (ix || iy) calc_edge[3] = *(PlaneZ(rear,iy)); \ else \ { LOGcalcNode(iz,vnode4,vnode0); calc_edge[3] = nvet; \ VX[nvet] = ix; VY[nvet] = iy; VZ[nvet] = cp; \ nvet++; \ } \ } \ if (exist_edge_4) /* cube vertex 1-3 */ \ { if (iz) calc_edge[4] = *(PlaneY(bellow,ix+1,iy)); \ else \ { calcNode(iy,vnode3,vnode1); calc_edge[4] = nvet; \ VX[nvet] = ix+1; VY[nvet] = cp; VZ[nvet] = iz; \ *(PlaneY(bellow,ix+1,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_5) /* cube vertex 1-5 */ \ { if (iy) calc_edge[5] = *(PlaneZ(front,iy)); \ else \ { LOGcalcNode(iz,vnode5,vnode1); calc_edge[5] = nvet; \ VX[nvet] = ix+1; VY[nvet] = iy; VZ[nvet] = cp; \ *(PlaneZ(front,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_6) /* cube vertex 2-3 */ \ { if (iz) calc_edge[6] = *(PlaneX(bellow,ix,iy+1)); \ else \ { calcNode(ix,vnode3,vnode2); calc_edge[6] = nvet; \ VX[nvet] = cp; VY[nvet] = iy+1; VZ[nvet] = iz; \ *(PlaneX(bellow,ix,iy+1)) = nvet; nvet++; \ } \ } \ if (exist_edge_7) /* cube vertex 2-6 */ \ { if (ix) calc_edge[7] = *(PlaneZ(rear,iy+1)); \ else \ { LOGcalcNode(iz,vnode6,vnode2); calc_edge[7] = nvet; \ VX[nvet] = ix; VY[nvet] = iy+1; VZ[nvet] = cp; \ *(PlaneZ(rear,iy+1)) = nvet; nvet++; \ } \ } \ if (exist_edge_8) /* cube vertex 3-7 */ \ { LOGcalcNode(iz,vnode7,vnode3); calc_edge[8] = nvet; \ VX[nvet] = ix+1; VY[nvet] = iy+1; VZ[nvet] = cp; \ *(PlaneZ(front,iy+1)) = nvet; nvet++; \ } \ if (exist_edge_9) /* cube vertex 4-5 */ \ { if (iy) calc_edge[9] = *(PlaneX(above,ix,iy)); \ else \ { calcNode(ix,vnode5,vnode4); calc_edge[9] = nvet; \ VX[nvet] = cp; VY[nvet] = iy; VZ[nvet] = iz+1; \ *(PlaneX(above,ix,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_A) /* cube vertex 4-6 */ \ { if (ix) calc_edge[10] = *(PlaneY(above,ix,iy)); \ else \ { calcNode(iy,vnode6,vnode4); calc_edge[10] = nvet; \ VX[nvet] = ix; VY[nvet] = cp; VZ[nvet] = iz+1; \ *(PlaneY(above,ix,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_B) /* cube vertex 5-7 */ \ { calcNode(iy,vnode7,vnode5); calc_edge[11] = nvet; \ VX[nvet] = ix+1; VY[nvet] = cp; VZ[nvet] = iz+1; \ *(PlaneY(above,ix+1,iy)) = nvet; nvet++; \ } \ if (exist_edge_C) /* cube vertex 6-7 */ \ { calcNode(ix,vnode7,vnode6); calc_edge[12] = nvet; \ VX[nvet] = cp; VY[nvet] = iy+1; VZ[nvet] = iz+1; \ *(PlaneX(above,ix,iy+1)) = nvet; nvet++; \ } \ } #endif #define INV_VAL INVALID_VALUE #define find_vertex_invalid_cube(nc) \ { ptFLAG[nc] &= 0x1FF; \ if (exist_polygon_in_cube(nc)) \ { if (exist_edge_4) /* cube vertex 1-3 */ \ { if (!iz && vnode3 < INV_VAL && vnode1 < INV_VAL) \ { calcNode(iy,vnode3,vnode1); \ VX[nvet] = ix+1; VY[nvet] = cp; VZ[nvet] = iz; \ *(PlaneY(bellow,ix+1,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_5) /* cube vertex 1-5 */ \ { if (!iy && vnode5 < INV_VAL && vnode1 < INV_VAL) \ { calcNode(iz,vnode5,vnode1); \ VX[nvet] = ix+1; VY[nvet] = iy; VZ[nvet] = cp; \ *(PlaneZ(front,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_6) /* cube vertex 2-3 */ \ { if (!iz && vnode3 < INV_VAL && vnode2 < INV_VAL) \ { calcNode(ix,vnode3,vnode2); \ VX[nvet] = cp; VY[nvet] = iy+1; VZ[nvet] = iz; \ *(PlaneX(bellow,ix,iy+1)) = nvet; nvet++; \ } \ } \ if (exist_edge_7) /* cube vertex 2-6 */ \ { if (!ix && vnode6 < INV_VAL && vnode2 < INV_VAL) \ { calcNode(iz,vnode6,vnode2); \ VX[nvet] = ix; VY[nvet] = iy+1; VZ[nvet] = cp; \ *(PlaneZ(rear,iy+1)) = nvet; nvet++; \ } \ } \ if (exist_edge_8) /* cube vertex 3-7 */ \ { if (vnode7 < INV_VAL && vnode3 < INV_VAL) \ { calcNode(iz,vnode7,vnode3); \ VX[nvet] = ix+1; VY[nvet] = iy+1; VZ[nvet] = cp; \ *(PlaneZ(front,iy+1)) = nvet; nvet++; \ } \ } \ if (exist_edge_9) /* cube vertex 4-5 */ \ { if (!iy && vnode5 < INV_VAL && vnode4 < INV_VAL) \ { calcNode(ix,vnode5,vnode4); \ VX[nvet] = cp; VY[nvet] = iy; VZ[nvet] = iz+1; \ *(PlaneX(above,ix,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_A) /* cube vertex 4-6 */ \ { if (!ix && vnode6 < INV_VAL && vnode4 < INV_VAL) \ { calcNode(iy,vnode6,vnode4); \ VX[nvet] = ix; VY[nvet] = cp; VZ[nvet] = iz+1; \ *(PlaneY(above,ix,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_B) /* cube vertex 5-7 */ \ { if (vnode7 < INV_VAL && vnode5 < INV_VAL) \ { calcNode(iy,vnode7,vnode5); \ VX[nvet] = ix+1; VY[nvet] = cp; VZ[nvet] = iz+1;\ *(PlaneY(above,ix+1,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_C) /* cube vertex 6-7 */ \ { if (vnode7 < INV_VAL && vnode6 < INV_VAL) \ { calcNode(ix,vnode7,vnode6); \ VX[nvet] = cp; VY[nvet] = iy+1; VZ[nvet] = iz+1;\ *(PlaneX(above,ix,iy+1)) = nvet; nvet++; \ } \ } \ } \ } #ifdef LEVELTYPES /* This macro is the same as find_vertex_invalid_cube, except for interpolation in the Z-coordinate. This is performed logaritmic by calling the macro LOGcalcNode. */ #define LOGfind_vertex_invalid_cube(nc) \ { ptFLAG[nc] &= 0x1FF; \ if (exist_polygon_in_cube(nc)) \ { if (exist_edge_4) /* cube vertex 1-3 */ \ { if (!iz && vnode3 < INV_VAL && vnode1 < INV_VAL) \ { calcNode(iy,vnode3,vnode1); \ VX[nvet] = ix+1; VY[nvet] = cp; VZ[nvet] = iz; \ *(PlaneY(bellow,ix+1,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_5) /* cube vertex 1-5 */ \ { if (!iy && vnode5 < INV_VAL && vnode1 < INV_VAL) \ { LOGcalcNode(iz,vnode5,vnode1); \ VX[nvet] = ix+1; VY[nvet] = iy; VZ[nvet] = cp; \ *(PlaneZ(front,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_6) /* cube vertex 2-3 */ \ { if (!iz && vnode3 < INV_VAL && vnode2 < INV_VAL) \ { calcNode(ix,vnode3,vnode2); \ VX[nvet] = cp; VY[nvet] = iy+1; VZ[nvet] = iz; \ *(PlaneX(bellow,ix,iy+1)) = nvet; nvet++; \ } \ } \ if (exist_edge_7) /* cube vertex 2-6 */ \ { if (!ix && vnode6 < INV_VAL && vnode2 < INV_VAL) \ { LOGcalcNode(iz,vnode6,vnode2); \ VX[nvet] = ix; VY[nvet] = iy+1; VZ[nvet] = cp; \ *(PlaneZ(rear,iy+1)) = nvet; nvet++; \ } \ } \ if (exist_edge_8) /* cube vertex 3-7 */ \ { if (vnode7 < INV_VAL && vnode3 < INV_VAL) \ { LOGcalcNode(iz,vnode7,vnode3); \ VX[nvet] = ix+1; VY[nvet] = iy+1; VZ[nvet] = cp; \ *(PlaneZ(front,iy+1)) = nvet; nvet++; \ } \ } \ if (exist_edge_9) /* cube vertex 4-5 */ \ { if (!iy && vnode5 < INV_VAL && vnode4 < INV_VAL) \ { calcNode(ix,vnode5,vnode4); \ VX[nvet] = cp; VY[nvet] = iy; VZ[nvet] = iz+1; \ *(PlaneX(above,ix,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_A) /* cube vertex 4-6 */ \ { if (!ix && vnode6 < INV_VAL && vnode4 < INV_VAL) \ { calcNode(iy,vnode6,vnode4); \ VX[nvet] = ix; VY[nvet] = cp; VZ[nvet] = iz+1; \ *(PlaneY(above,ix,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_B) /* cube vertex 5-7 */ \ { if (vnode7 < INV_VAL && vnode5 < INV_VAL) \ { calcNode(iy,vnode7,vnode5); \ VX[nvet] = ix+1; VY[nvet] = cp; VZ[nvet] = iz+1;\ *(PlaneY(above,ix+1,iy)) = nvet; nvet++; \ } \ } \ if (exist_edge_C) /* cube vertex 6-7 */ \ { if (vnode7 < INV_VAL && vnode6 < INV_VAL) \ { calcNode(ix,vnode7,vnode6); \ VX[nvet] = cp; VY[nvet] = iy+1; VZ[nvet] = iz+1;\ *(PlaneX(above,ix,iy+1)) = nvet; nvet++; \ } \ } \ } \ } #endif #define num_pol_in_cube(nc) pol_edges[ptFLAG[nc]][1] #define num_edges_per_pol(nc) pol_edges[ptFLAG[nc]][2] #define vet_edges_pol(nc) pol_edges[ptFLAG[nc]][4] #define fill_Vert_f_Pol( nc ) \ { kk = num_edges_per_pol(nc); \ ptl = &vet_edges_pol(nc); \ pa = pvp; \ for (ii=0; ii>= 4; pa += 7; \ } \ } #define update_data_structure( nc ) \ { kk = num_edges_per_pol(nc); \ nn = num_pol_in_cube(nc); \ for (ii=0; ii>= 4; pvp += 7; cpl++; \ } \ } #define update_data_structure_with_print( nc ) \ { kk = num_edges_per_pol(nc); \ nn = num_pol_in_cube(nc); \ for (ii=0; ii>= 4; pvp += 7; cpl++; \ } \ } #define swap_planes(a,p1,p2) \ { caseA = p1; p1 = p2; p2 = caseA; } static int marching( Context ctx, float ptGRID[], int xdim, int ydim, int zdim, int ptFLAG[], float VX[], float VY[], float VZ[], int NVERTICE, double isovalue, int npolygons, int Pol_f_Vert[], int Vert_f_Pol[] ) { int *ixPlane[2], *iyPlane[2], *izPlane[2]; /* Plane Address Parameters */ int ix, iy, iz, caseA, above, bellow, front, rear, mm, nn; int *ptl, ii, jj, kk, ncube, cpl, pvp, pa, ve, calc_edge[13]; float cp, *pt; double nodeDiff; int xdim_x_ydim = xdim*ydim; int nvertex; bellow = rear = 0; above = front = 1; /* Initialize the Auxiliar Arrays of Pointers */ ix = 9 * (npolygons*2 + 50); iy = 7 * npolygons; ii = ix + iy; /*$dir vector */ for (jj=0; jj NVERTICE) goto end; if (valid_cube(ncube)) { fill_Vert_f_Pol(ncube); #ifdef LEVELTYPES switch(LevelType) { case PRESSURE_LEVELS: LOGfind_vertex(); break; default: find_vertex(); break; } #else find_vertex(); #endif update_data_structure_with_print(ncube); } else { #ifdef LEVELTYPES switch(LevelType) { case PRESSURE_LEVELS: LOGfind_vertex_invalid_cube(ncube); break; default: find_vertex_invalid_cube(ncube); break; } #else find_vertex_invalid_cube(ncube); #endif } } else print_no_polygon(); ncube++; pt++; } swap_planes(Z,rear,front); pt++; } swap_planes(XY,bellow,above); pt += ydim; } #else for ( iz = 0; iz < zdim - 1; iz++ ) { for ( ix = 0; ix < xdim - 1; ix++ ) { for ( iy = 0; iy < ydim - 1; iy++ ) { if (exist_polygon_in_cube(ncube)) { if ((nvertex + 12) > NVERTICE) goto end; if (valid_cube(ncube)) { fill_Vert_f_Pol(ncube); #ifdef LEVELTYPES switch(LevelType) { case PRESSURE_LEVELS: LOGfind_vertex(); break; default: find_vertex(); break; } #else find_vertex(); #endif update_data_structure(ncube); } else { #ifdef LEVELTYPES switch(LevelType) { case PRESSURE_LEVELS: LOGfind_vertex_invalid_cube(ncube); break; default: find_vertex_invalid_cube(ncube); break; } #else find_vertex_invalid_cube(ncube); #endif } } ncube++; pt++; } swap_planes(Z,rear,front); pt++; } swap_planes(XY,bellow,above); pt += ydim; } #endif end: xfree ( ctx, (char *) ixPlane[0]); return nvertex; } /***********************************************************/ /**** NORMALS ****/ /***********************************************************/ #define GV(x,y,z) ptGRID[(int)y + (int)x*ydim + (int)z*xdim_x_ydim] #define Vertex_k_Pol_i pVP[k] #define V_f_P Vert_f_Pol #define adjust_normal_by_gradiente(v,x,y,z) \ { if (VX[v[0]]==VX[v[1]] || VX[v[0]]==VX[v[2]]) \ ix = (int)VX[v[0]]; \ else \ ix = (int)VX[v[1]]; \ if (VY[v[0]]==VY[v[1]] || VY[v[0]]==VY[v[2]]) \ iy = (int)VY[v[0]]; \ else \ iy = (int)VY[v[1]]; \ if (VZ[v[0]]==VZ[v[1]] || VZ[v[0]]==VZ[v[2]]) \ iz = (int)VZ[v[0]]; \ else \ iz = (int)VZ[v[1]]; \ i1 = ix; \ if (i1 != ixb) i2 = i1 + 1; \ else { i2 = i1; i1--; } \ if (GV(i2,iy,iz) >= INVALID_VALUE) { i2--; i1--; } \ x = GV(i2,iy,iz) - GV(i1,iy,iz); \ i1 = iy; \ if (i1 != iyb) i2 = i1 + 1; \ else { i2 = i1; i1--; } \ if (GV(ix,i2,iz) >= INVALID_VALUE) { i2--; i1--; } \ y = GV(ix,i2,iz) - GV(ix,i1,iz); \ i1 = iz; \ if (i1 != izb) i2 = i1 + 1; \ else { i2 = i1; i1--; } \ if (GV(ix,iy,i2) >= INVALID_VALUE) { i2--; i1--; } \ z = GV(ix,iy,i2) - GV(ix,iy,i1); \ a = (x*x + y*y + z*z); \ if (a > 0.) { x /= a; y /= a; z /= a; } \ } static void normals( float ptGRID[], int xdim, int ydim, int zdim, float VX[], float VY[], float VZ[], float NX[], float NY[], float NZ[], int nvertex, int npolygons, float Pnx[], float Pny[], float Pnz[], float NxA[], float NxB[], float NyA[], float NyB[], float NzA[], float NzB[], double arX, double arY, double arZ, int Pol_f_Vert[], int Vert_f_Pol[] ) { int i, k, iv[3], n; int i1, i2, ix, iy, iz, ixb, iyb, izb, *pVP; int max_vert_per_pol, swap_flag; float x, y, z, a, minimum_area, len; int xdim_x_ydim = xdim*ydim; ixb = xdim-1; iyb = ydim-1; izb = zdim-1; minimum_area = max(1.e-4,EPS_0); /* Calculate maximum number of vertices per polygon */ k = 6; n = 7*npolygons; while ( TRUE ) { for (i=k+7; i Vert_f_Pol[k]) break; if (i >= n) break; k = i; } max_vert_per_pol = Vert_f_Pol[k]; /* Calculate the Normals vector components for each Polygon */ pVP = Vert_f_Pol; /*$dir vector */ for ( i=0; i0) { /* check for valid polygon added by BEP 2-13-92 */ NxA[i] = VX[pVP[1]] - VX[pVP[0]]; NyA[i] = VY[pVP[1]] - VY[pVP[0]]; NzA[i] = VZ[pVP[1]] - VZ[pVP[0]]; } pVP += 7; } swap_flag = 0; for ( k = 2; k < max_vert_per_pol; k++ ) { pVP = Vert_f_Pol; if (swap_flag==0) { /*$dir no_recurrence */ /* Vectorized */ for ( i=0; i= 0) { NxB[i] = VX[pVP[k]] - VX[pVP[0]]; NyB[i] = VY[pVP[k]] - VY[pVP[0]]; NzB[i] = VZ[pVP[k]] - VZ[pVP[0]]; Pnx[i] = NyA[i]*NzB[i] - NzA[i]*NyB[i]; Pny[i] = NzA[i]*NxB[i] - NxA[i]*NzB[i]; Pnz[i] = NxA[i]*NyB[i] - NyA[i]*NxB[i]; NxA[i] = Pnx[i]*Pnx[i] + Pny[i]*Pny[i] + Pnz[i]*Pnz[i]; if (NxA[i] > minimum_area) { Pnx[i] /= NxA[i]; Pny[i] /= NxA[i]; Pnz[i] /= NxA[i]; } } pVP += 7; } } else { /* swap_flag!=0 */ /*$dir no_recurrence */ /* Vectorized */ for ( i=0; i= 0) { NxA[i] = VX[pVP[k]] - VX[pVP[0]]; NyA[i] = VY[pVP[k]] - VY[pVP[0]]; NzA[i] = VZ[pVP[k]] - VZ[pVP[0]]; Pnx[i] = NyB[i]*NzA[i] - NzB[i]*NyA[i]; Pny[i] = NzB[i]*NxA[i] - NxB[i]*NzA[i]; Pnz[i] = NxB[i]*NyA[i] - NyB[i]*NxA[i]; NxB[i] = Pnx[i]*Pnx[i] + Pny[i]*Pny[i] + Pnz[i]*Pnz[i]; if (NxB[i] > minimum_area) { Pnx[i] /= NxB[i]; Pny[i] /= NxB[i]; Pnz[i] /= NxB[i]; } } pVP += 7; } } pVP = Vert_f_Pol; /* This Loop be Vectorized */ for ( i=0; i= 0) { iv[0] = pVP[0]; iv[1] = pVP[k-1]; iv[2] = pVP[k]; if (NxA[i] > minimum_area) { x = Pnx[i]; y = Pny[i]; z = Pnz[i]; } else adjust_normal_by_gradiente(iv,x,y,z); /* Update the origin vertex */ NX[iv[0]] += x; NY[iv[0]] += y; NZ[iv[0]] += z; /* Update the vertex that defines the first vector */ NX[iv[1]] += x; NY[iv[1]] += y; NZ[iv[1]] += z; /* Update the vertex that defines the second vector */ NX[iv[2]] += x; NY[iv[2]] += y; NZ[iv[2]] += z; } pVP += 7; } swap_flag = (swap_flag ? 0 : 1 ); } /* Apply Aspect Ratio in the Normals */ if (arX != 1.0) for (i=0; i EPS_0) { NX[i] /= len; NY[i] /= len; NZ[i] /= len; } } } /**********************************************************/ /**** TRIANGLES ****/ /**********************************************************/ #define P_f_V Pol_f_Vert #define V_f_P Vert_f_Pol #define swap(a,b) { i = a; a = b; b = i; } #define find_unselected_pol(p) \ { for (p=last_pol; p=0 && vb>=0) { /* check added by BEP 2-13-92 */ \ i=va*9; k=i+P_f_V[i+8]; j=vb*9; m=j+P_f_V[j+8]; \ while (i>0 && j>0 && i= 0) break; } /* insert polygon alone */ if (npol < 0) { ptT = NTAB + STAB[Nvt-3]; if (iST > 0) { Tri_Stripe[iST] = Tri_Stripe[iST-1]; iST++; Tri_Stripe[iST++] = Vt[*ptT]; } else f_line_conection = 1; /* WLH 3-9-95 added */ for (ii=0; ii EPS_0) for (i=0; i EPS_0) for (i=0; i EPS_0) for (i=0; i number of ocorrunces of the "flag case f" - nvertex -> number of vertex generated - npolygons -> number of polygons generated - num_cubes -> number of cubes in the matrix */ #ifdef DEBUG #define TOTFLAGS 268 static void statistics( double isovalue, int num_cubes ) { int i, total, vt_flags[TOTFLAGS]; float nc; fprintf( output, "\n== Statistics ==\n" ); fprintf( output, "Isovalue: %9.3e\n", isovalue ); fprintf( output, "Nodes generated: %7d\n", nvertex ); fprintf( output, "Polygons generated: %7d\n", npolygons ); fprintf( output, "Cubes in matrix: %7d\n", num_cubes ); nc = num_cubes / 100.0; for (i=0; i %02x): %7d (%5.2f%%)\n", i, i, pol_edges[i][0], vt_flags[i], (float)vt_flags[i]/nc ); } fprintf( output, " total : %7d (100.00%%)\n\n", num_cubes); fprintf( output, " total cases #0: %7d (%5.2f%%)\n", total, (float)total / nc ); return; } /* ^ ^ ----- end of the statistics routine */ #endif