BICPL: marching_tetra.c Source File

00001 00002 #include <volume_io/internal_volume_io.h> 00003 #include <bicpl/marching.h> 00004 00005 #define MAX_POLYGONS_PER_VOXEL 5 00006 #define MAX_INDICES_PER_VOXEL 20 00007 00008 typedef struct 00009 { 00010 int n_polygons; 00011 int *sizes; 00012 voxel_point_type *indices; 00013 } case_struct; 00014 00015 typedef enum { PLUS_FLAG, MINUS_FLAG, MAX_CASES } 00016 Case_types; 00017 00018 private case_struct cases[2][2][2][2][2][2][2][2][2][2][2]; 00019 00020 private BOOLEAN initialized = FALSE; 00021 00022 /* ---------------- private prototypes ------------------------- */ 00023 00024 private void create_marching_cubes_lookup( void ); 00025 private void create_case( 00026 int x, 00027 int y, 00028 int z, 00029 Case_types case_flags[2][2][2], 00030 case_struct *case_info ); 00031 00032 /* ------------------------------------------------------------- */ 00033 00034 private void check_initialized( void ) 00035 { 00036 if( !initialized ) 00037 { 00038 create_marching_cubes_lookup(); 00039 00040 initialized = TRUE; 00041 } 00042 } 00043 00044 private void create_marching_cubes_lookup( void ) 00045 { 00046 int x, y, z; 00047 Case_types case_flags[2][2][2]; 00048 00049 for_less( x, 0, 2 ) 00050 for_less( y, 0, 2 ) 00051 for_less( z, 0, 2 ) 00052 for_enum( case_flags[0][0][0], MAX_CASES, Case_types ) 00053 for_enum( case_flags[0][0][1], MAX_CASES, Case_types ) 00054 for_enum( case_flags[0][1][0], MAX_CASES, Case_types ) 00055 for_enum( case_flags[0][1][1], MAX_CASES, Case_types ) 00056 for_enum( case_flags[1][0][0], MAX_CASES, Case_types ) 00057 for_enum( case_flags[1][0][1], MAX_CASES, Case_types ) 00058 for_enum( case_flags[1][1][0], MAX_CASES, Case_types ) 00059 for_enum( case_flags[1][1][1], MAX_CASES, Case_types ) 00060 { 00061 create_case( x, y, z, case_flags, 00062 &cases[x][y][z][case_flags[0][0][0]] 00063 [case_flags[0][0][1]] 00064 [case_flags[0][1][0]] 00065 [case_flags[0][1][1]] 00066 [case_flags[1][0][0]] 00067 [case_flags[1][0][1]] 00068 [case_flags[1][1][0]] 00069 [case_flags[1][1][1]] ); 00070 } 00071 } 00072 00073 private int get_tetra_polygon( 00074 Case_types tetra_flags[], 00075 int tetra_indices[][2] ) 00076 { 00077 static int sizes[5] = { 0, 3, 4, 3, 0 }; 00078 static int indices[5][4][2] = { 00079 { {0,0} }, 00080 { {0,1}, {0,2}, {0,3} }, 00081 { {0,2}, {0,3}, {1,3}, {1,2} }, 00082 { {0,3}, {1,3}, {2,3} }, 00083 { {0,0} } 00084 }; 00085 int i, n_minus; 00086 00087 n_minus = 0; 00088 for_less( i, 0, 4 ) 00089 { 00090 if( tetra_flags[i] == MINUS_FLAG ) 00091 ++n_minus; 00092 } 00093 00094 for_less( i, 0, sizes[n_minus] ) 00095 { 00096 tetra_indices[i][0] = indices[n_minus][i][0]; 00097 tetra_indices[i][1] = indices[n_minus][i][1]; 00098 } 00099 00100 return( sizes[n_minus] ); 00101 } 00102 00103 private int get_tetra_case( 00104 Case_types case_flags[2][2][2], 00105 int indices[][3], 00106 int edge_indices[][2][3] ) 00107 { 00108 int ord, rank, best_ord, best_rank, vertex; 00109 int tetra_indices[4][2], x, y, z, p, e, size; 00110 Case_types tetra_flags[4]; 00111 static int orders[12][4] = { 00112 { 0, 1, 2, 3 }, 00113 { 0, 2, 3, 1 }, 00114 { 0, 3, 1, 2 }, 00115 { 1, 2, 0, 3 }, 00116 { 1, 0, 3, 2 }, 00117 { 1, 3, 2, 0 }, 00118 { 2, 0, 1, 3 }, 00119 { 2, 1, 3, 0 }, 00120 { 2, 3, 0, 1 }, 00121 { 3, 0, 2, 1 }, 00122 { 3, 2, 1, 0 }, 00123 { 3, 1, 0, 2 } 00124 }; 00125 00126 best_rank = 0; 00127 00128 for_less( ord, 0, 12 ) 00129 { 00130 rank = 0; 00131 for_less( vertex, 0, 4 ) 00132 { 00133 x = indices[orders[ord][vertex]][0]; 00134 y = indices[orders[ord][vertex]][1]; 00135 z = indices[orders[ord][vertex]][2]; 00136 if( case_flags[x][y][z] == PLUS_FLAG ) 00137 rank |= (1 << (3 - vertex)); 00138 } 00139 00140 if( ord == 0 || rank < best_rank ) 00141 { 00142 best_rank = rank; 00143 best_ord = ord; 00144 } 00145 } 00146 00147 for_less( vertex, 0, 4 ) 00148 { 00149 x = indices[orders[best_ord][vertex]][0]; 00150 y = indices[orders[best_ord][vertex]][1]; 00151 z = indices[orders[best_ord][vertex]][2]; 00152 tetra_flags[vertex] = case_flags[x][y][z]; 00153 } 00154 00155 size = get_tetra_polygon( tetra_flags, tetra_indices ); 00156 00157 for_less( p, 0, size ) 00158 { 00159 for_less( e, 0, 2 ) 00160 { 00161 vertex = orders[best_ord][tetra_indices[p][e]]; 00162 00163 edge_indices[p][e][0] = indices[vertex][0]; 00164 edge_indices[p][e][1] = indices[vertex][1]; 00165 edge_indices[p][e][2] = indices[vertex][2]; 00166 } 00167 } 00168 00169 return( size ); 00170 } 00171 00172 private int lookup_case( 00173 Case_types case_flags[2][2][2], 00174 int sizes[], 00175 int edge_indices[][2][3] ) 00176 { 00177 int n_polygons, ind, size, tetra; 00178 static int indices[5][4][3] = { 00179 { { 0, 0, 0 }, 00180 { 0, 0, 1 }, 00181 { 0, 1, 1 }, 00182 { 1, 0, 1 } }, 00183 { { 0, 0, 0 }, 00184 { 1, 0, 0 }, 00185 { 1, 0, 1 }, 00186 { 1, 1, 0 } }, 00187 { { 0, 1, 1 }, 00188 { 1, 1, 1 }, 00189 { 1, 1, 0 }, 00190 { 1, 0, 1 } }, 00191 { { 0, 1, 0 }, 00192 { 0, 1, 1 }, 00193 { 1, 1, 0 }, 00194 { 0, 0, 0 } }, 00195 { { 0, 0, 0 }, 00196 { 0, 1, 1 }, 00197 { 1, 1, 0 }, 00198 { 1, 0, 1 } } 00199 }; 00200 00201 n_polygons = 0; 00202 ind = 0; 00203 for_less( tetra, 0, 5 ) 00204 { 00205 size = get_tetra_case( case_flags, indices[tetra], &edge_indices[ind] ); 00206 if( size > 0 ) 00207 { 00208 sizes[n_polygons] = size; 00209 ++n_polygons; 00210 ind += size; 00211 } 00212 } 00213 00214 return( n_polygons ); 00215 } 00216 00217 static int offsets[N_MARCHING_TETRA_EDGES][3] = { 00218 { 1, 0, 0 }, 00219 { 0, 1, 0 }, 00220 { 0, 0, 1 }, 00221 { 0, 1, 1 }, 00222 { 1, 0, 1 }, 00223 { 1, 1, 0 }, 00224 { 1, 1, 1 }, 00225 { 1, 0, -1 }, 00226 { 0, 1, -1 }, 00227 { 1, -1, 0 }, 00228 { 0, -1, -1 } 00229 }; 00230 00231 public void translate_to_edge_index( 00232 int x1, 00233 int y1, 00234 int z1, 00235 int x2, 00236 int y2, 00237 int z2, 00238 voxel_point_type *edge_point ) 00239 { 00240 int tmp, edge; 00241 00242 if( x2 < x1 || (x1 == x2 && y2 < y1) || (x1 == x2 && y1 == y2 && z2 < z1) ) 00243 { 00244 tmp = x1; 00245 x1 = x2; 00246 x2 = tmp; 00247 tmp = y1; 00248 y1 = y2; 00249 y2 = tmp; 00250 tmp = z1; 00251 z1 = z2; 00252 z2 = tmp; 00253 } 00254 00255 for_less( edge, 0, SIZEOF_STATIC_ARRAY(offsets) ) 00256 { 00257 if( x1 + offsets[edge][X] == x2 && 00258 y1 + offsets[edge][Y] == y2 && 00259 z1 + offsets[edge][Z] == z2 ) 00260 break; 00261 } 00262 00263 if( edge >= SIZEOF_STATIC_ARRAY(offsets) ) 00264 handle_internal_error( "edge_intersected" ); 00265 00266 edge_point->edge_intersected = edge; 00267 edge_point->coord[0] = x1; 00268 edge_point->coord[1] = y1; 00269 edge_point->coord[2] = z1; 00270 } 00271 00272 public void translate_from_edge_index( 00273 int edge_index, 00274 int offset[] ) 00275 { 00276 offset[0] = offsets[edge_index][0]; 00277 offset[1] = offsets[edge_index][1]; 00278 offset[2] = offsets[edge_index][2]; 00279 } 00280 00281 private void create_case( 00282 int x, 00283 int y, 00284 int z, 00285 Case_types case_flags[2][2][2], 00286 case_struct *case_info ) 00287 { 00288 int i, j, k, tx, ty, tz, x1, y1, z1, x2, y2, z2; 00289 int poly, n_indices, ind, vertex; 00290 voxel_point_type tmp; 00291 int translation_indices[2][2][2][3]; 00292 int edge_indices[MAX_INDICES_PER_VOXEL][2][3]; 00293 Case_types transformed_case[2][2][2]; 00294 int sizes[MAX_POLYGONS_PER_VOXEL]; 00295 BOOLEAN left_handed; 00296 00297 for_less( i, 0, 2 ) 00298 for_less( j, 0, 2 ) 00299 for_less( k, 0, 2 ) 00300 { 00301 tx = i ^ x; 00302 ty = j ^ y; 00303 tz = k ^ z; 00304 translation_indices[i][j][k][0] = tx; 00305 translation_indices[i][j][k][1] = ty; 00306 translation_indices[i][j][k][2] = tz; 00307 transformed_case[i][j][k] = case_flags[tx][ty][tz]; 00308 } 00309 00310 case_info->n_polygons = lookup_case( transformed_case, 00311 sizes, edge_indices ); 00312 00313 if( case_info->n_polygons == 0 ) 00314 return; 00315 00316 ALLOC( case_info->sizes, case_info->n_polygons ); 00317 n_indices = 0; 00318 for_less( poly, 0, case_info->n_polygons ) 00319 { 00320 case_info->sizes[poly] = sizes[poly]; 00321 n_indices += sizes[poly]; 00322 } 00323 00324 ALLOC( case_info->indices, n_indices ); 00325 00326 for_less( ind, 0, n_indices ) 00327 { 00328 tx = edge_indices[ind][0][X]; 00329 ty = edge_indices[ind][0][Y]; 00330 tz = edge_indices[ind][0][Z]; 00331 x1 = translation_indices[tx][ty][tz][X]; 00332 y1 = translation_indices[tx][ty][tz][Y]; 00333 z1 = translation_indices[tx][ty][tz][Z]; 00334 00335 tx = edge_indices[ind][1][X]; 00336 ty = edge_indices[ind][1][Y]; 00337 tz = edge_indices[ind][1][Z]; 00338 x2 = translation_indices[tx][ty][tz][X]; 00339 y2 = translation_indices[tx][ty][tz][Y]; 00340 z2 = translation_indices[tx][ty][tz][Z]; 00341 00342 translate_to_edge_index( x1, y1, z1, x2, y2, z2, 00343 &case_info->indices[ind] ); 00344 } 00345 00346 left_handed = ((x + y + z) % 2) == 1; 00347 00348 if( left_handed ) 00349 { 00350 ind = 0; 00351 for_less( poly, 0, case_info->n_polygons ) 00352 { 00353 for_less( vertex, 0, sizes[poly]/2 ) 00354 { 00355 tmp = case_info->indices[ind+vertex]; 00356 case_info->indices[ind+vertex] = 00357 case_info->indices[ind+sizes[poly]-1-vertex]; 00358 case_info->indices[ind+sizes[poly]-1-vertex] = tmp; 00359 } 00360 ind += sizes[poly]; 00361 } 00362 } 00363 } 00364 00365 private void delete_case( 00366 case_struct *case_info ); 00367 00368 public int get_tetra_isosurface_polygons( 00369 int x, 00370 int y, 00371 int z, 00372 Real corners[2][2][2], 00373 Real isovalue, 00374 int *sizes[], 00375 voxel_point_type *points[] ) 00376 { 00377 int c0, c1, c2, c3, c4, c5, c6, c7, xc, yc, zc; 00378 case_struct *voxel_case; 00379 00380 check_initialized(); 00381 00382 c0 = (corners[0][0][0] <= isovalue); 00383 c1 = (corners[0][0][1] <= isovalue); 00384 c2 = (corners[0][1][0] <= isovalue); 00385 c3 = (corners[0][1][1] <= isovalue); 00386 c4 = (corners[1][0][0] <= isovalue); 00387 c5 = (corners[1][0][1] <= isovalue); 00388 c6 = (corners[1][1][0] <= isovalue); 00389 c7 = (corners[1][1][1] <= isovalue); 00390 00391 xc = x & 1; 00392 yc = y & 1; 00393 zc = z & 1; 00394 00395 voxel_case = &cases[xc][yc][zc][c0][c1][c2][c3][c4][c5][c6][c7]; 00396 00397 #ifdef DEBUG 00398 { 00399 Case_types case_flags[2][2][2]; 00400 00401 case_flags[0][0][0] = (Case_types) c0; 00402 case_flags[0][0][1] = (Case_types) c1; 00403 case_flags[0][1][0] = (Case_types) c2; 00404 case_flags[0][1][1] = (Case_types) c3; 00405 case_flags[1][0][0] = (Case_types) c4; 00406 case_flags[1][0][1] = (Case_types) c5; 00407 case_flags[1][1][0] = (Case_types) c6; 00408 case_flags[1][1][1] = (Case_types) c7; 00409 delete_case( voxel_case ); 00410 create_case( xc, yc, zc, case_flags, voxel_case ); 00411 } 00412 #endif 00413 00414 *sizes = voxel_case->sizes; 00415 *points = voxel_case->indices; 00416 00417 return( voxel_case->n_polygons ); 00418 } 00419 00420 private void delete_case( 00421 case_struct *case_info ) 00422 { 00423 if( case_info->n_polygons > 0 ) 00424 { 00425 FREE( case_info->sizes ); 00426 FREE( case_info->indices ); 00427 } 00428 } 00429 00430 public void delete_tetra_marching_cubes_table( void ) 00431 { 00432 int x, y, z; 00433 Case_types c0, c1, c2, c3, c4, c5, c6, c7; 00434 00435 if( initialized ) 00436 { 00437 for_less( x, 0, 2 ) 00438 for_less( y, 0, 2 ) 00439 for_less( z, 0, 2 ) 00440 for_enum( c0, MAX_CASES, Case_types ) 00441 for_enum( c1, MAX_CASES, Case_types ) 00442 for_enum( c2, MAX_CASES, Case_types ) 00443 for_enum( c3, MAX_CASES, Case_types ) 00444 for_enum( c4, MAX_CASES, Case_types ) 00445 for_enum( c5, MAX_CASES, Case_types ) 00446 for_enum( c6, MAX_CASES, Case_types ) 00447 for_enum( c7, MAX_CASES, Case_types ) 00448 { 00449 delete_case( &cases[x][y][z][c0][c1][c2][c3][c4][c5][c6][c7] ); 00450 } 00451 } 00452 }