GRASS Programmer's Manual
6.4.2(2012)
|
00001 00022 #include <float.h> 00023 00024 #include <grass/gis.h> 00025 #include <grass/gstypes.h> 00026 00027 #include "mc33_table.h" 00028 00029 unsigned char m_case, m_config, m_subconfig; 00030 00039 int mc33_test_face(char face, float *v) 00040 { 00041 float A, B, C, D; 00042 00043 switch (face) { 00044 case -1: 00045 case 1: 00046 A = v[0]; 00047 B = v[4]; 00048 C = v[5]; 00049 D = v[1]; 00050 break; 00051 00052 case -2: 00053 case 2: 00054 A = v[1]; 00055 B = v[5]; 00056 C = v[6]; 00057 D = v[2]; 00058 break; 00059 00060 case -3: 00061 case 3: 00062 A = v[2]; 00063 B = v[6]; 00064 C = v[7]; 00065 D = v[3]; 00066 break; 00067 00068 case -4: 00069 case 4: 00070 A = v[3]; 00071 B = v[7]; 00072 C = v[4]; 00073 D = v[0]; 00074 break; 00075 00076 case -5: 00077 case 5: 00078 A = v[0]; 00079 B = v[3]; 00080 C = v[2]; 00081 D = v[1]; 00082 break; 00083 00084 case -6: 00085 case 6: 00086 A = v[4]; 00087 B = v[7]; 00088 C = v[6]; 00089 D = v[5]; 00090 break; 00091 00092 default: 00093 fprintf(stderr, "Invalid face code %d\n", face); 00094 A = B = C = D = 0; 00095 }; 00096 00097 return face * A * (A * C - B * D) >= 0; 00098 } 00099 00108 int mc33_test_interior(char s, float *v) 00109 { 00110 float t, At = 0, Bt = 0, Ct = 0, Dt = 0, a, b; 00111 char test = 0; 00112 char edge = -1; 00113 00114 switch (m_case) { 00115 case 4: 00116 case 10: 00117 a = (v[4] - v[0]) * (v[6] - v[2]) - (v[7] - v[3]) * (v[5] - v[1]); 00118 b = v[2] * (v[4] - v[0]) + v[0] * (v[6] - v[2]) 00119 - v[1] * (v[7] - v[3]) - v[3] * (v[5] - v[1]); 00120 t = -b / (2 * a); 00121 00122 if (t < 0 || t > 1) 00123 return s > 0; 00124 00125 At = v[0] + (v[4] - v[0]) * t; 00126 Bt = v[3] + (v[7] - v[3]) * t; 00127 Ct = v[2] + (v[6] - v[2]) * t; 00128 Dt = v[1] + (v[5] - v[1]) * t; 00129 break; 00130 00131 case 6: 00132 case 7: 00133 case 12: 00134 case 13: 00135 switch (m_case) { 00136 case 6: 00137 edge = cell_table[OFFSET_T6_1_1 + m_config].polys[0]; 00138 break; 00139 case 7: 00140 edge = cell_table[OFFSET_T7_4_1 + m_config].polys[13]; 00141 break; 00142 case 12: 00143 edge = cell_table[OFFSET_T12_2_S1 + m_config].polys[14]; 00144 break; 00145 case 13: 00146 edge = 00147 cell_table[OFFSET_T13_5_1 + m_subconfig + 00148 m_config * 4].polys[2]; 00149 break; 00150 } 00151 00152 switch (edge) { 00153 case 0: 00154 t = v[0] / (v[0] - v[1]); 00155 At = 0; 00156 Bt = v[3] + (v[2] - v[3]) * t; 00157 Ct = v[7] + (v[6] - v[7]) * t; 00158 Dt = v[4] + (v[5] - v[4]) * t; 00159 break; 00160 case 1: 00161 t = v[1] / (v[1] - v[2]); 00162 At = 0; 00163 Bt = v[0] + (v[3] - v[0]) * t; 00164 Ct = v[4] + (v[7] - v[4]) * t; 00165 Dt = v[5] + (v[6] - v[5]) * t; 00166 break; 00167 case 2: 00168 t = v[2] / (v[2] - v[3]); 00169 At = 0; 00170 Bt = v[1] + (v[0] - v[1]) * t; 00171 Ct = v[5] + (v[4] - v[5]) * t; 00172 Dt = v[6] + (v[7] - v[6]) * t; 00173 break; 00174 case 3: 00175 t = v[3] / (v[3] - v[0]); 00176 At = 0; 00177 Bt = v[2] + (v[1] - v[2]) * t; 00178 Ct = v[6] + (v[5] - v[6]) * t; 00179 Dt = v[7] + (v[4] - v[7]) * t; 00180 break; 00181 case 4: 00182 t = v[4] / (v[4] - v[5]); 00183 At = 0; 00184 Bt = v[7] + (v[6] - v[7]) * t; 00185 Ct = v[3] + (v[2] - v[3]) * t; 00186 Dt = v[0] + (v[1] - v[0]) * t; 00187 break; 00188 case 5: 00189 t = v[5] / (v[5] - v[6]); 00190 At = 0; 00191 Bt = v[4] + (v[7] - v[4]) * t; 00192 Ct = v[0] + (v[3] - v[0]) * t; 00193 Dt = v[1] + (v[2] - v[1]) * t; 00194 break; 00195 case 6: 00196 t = v[6] / (v[6] - v[7]); 00197 At = 0; 00198 Bt = v[5] + (v[4] - v[5]) * t; 00199 Ct = v[1] + (v[0] - v[1]) * t; 00200 Dt = v[2] + (v[3] - v[2]) * t; 00201 break; 00202 case 7: 00203 t = v[7] / (v[7] - v[4]); 00204 At = 0; 00205 Bt = v[6] + (v[5] - v[6]) * t; 00206 Ct = v[2] + (v[1] - v[2]) * t; 00207 Dt = v[3] + (v[0] - v[3]) * t; 00208 break; 00209 case 8: 00210 t = v[0] / (v[0] - v[4]); 00211 At = 0; 00212 Bt = v[3] + (v[7] - v[3]) * t; 00213 Ct = v[2] + (v[6] - v[2]) * t; 00214 Dt = v[1] + (v[5] - v[1]) * t; 00215 break; 00216 case 9: 00217 t = v[1] / (v[1] - v[5]); 00218 At = 0; 00219 Bt = v[0] + (v[4] - v[0]) * t; 00220 Ct = v[3] + (v[7] - v[3]) * t; 00221 Dt = v[2] + (v[6] - v[2]) * t; 00222 break; 00223 case 10: 00224 t = v[2] / (v[2] - v[6]); 00225 At = 0; 00226 Bt = v[1] + (v[5] - v[1]) * t; 00227 Ct = v[0] + (v[4] - v[0]) * t; 00228 Dt = v[3] + (v[7] - v[3]) * t; 00229 break; 00230 case 11: 00231 t = v[3] / (v[3] - v[7]); 00232 At = 0; 00233 Bt = v[2] + (v[6] - v[2]) * t; 00234 Ct = v[1] + (v[5] - v[1]) * t; 00235 Dt = v[0] + (v[4] - v[0]) * t; 00236 break; 00237 default: 00238 fprintf(stderr, "Invalid edge %d\n", edge); 00239 break; 00240 } 00241 break; 00242 00243 default: 00244 fprintf(stderr, "Invalid ambiguous case %d\n", m_case); 00245 break; 00246 } 00247 00248 if (At >= 0) 00249 test++; 00250 if (Bt >= 0) 00251 test += 2; 00252 if (Ct >= 0) 00253 test += 4; 00254 if (Dt >= 0) 00255 test += 8; 00256 00257 switch (test) { 00258 case 0: 00259 return s > 0; 00260 case 1: 00261 return s > 0; 00262 case 2: 00263 return s > 0; 00264 case 3: 00265 return s > 0; 00266 case 4: 00267 return s > 0; 00268 case 5: 00269 if (At * Ct < Bt * Dt) 00270 return s > 0; 00271 break; 00272 case 6: 00273 return s > 0; 00274 case 7: 00275 return s < 0; 00276 case 8: 00277 return s > 0; 00278 case 9: 00279 return s > 0; 00280 case 10: 00281 if (At * Ct >= Bt * Dt) 00282 return s > 0; 00283 break; 00284 case 11: 00285 return s < 0; 00286 case 12: 00287 return s > 0; 00288 case 13: 00289 return s < 0; 00290 case 14: 00291 return s < 0; 00292 case 15: 00293 return s < 0; 00294 } 00295 00296 return s < 0; 00297 } 00298 00307 int mc33_process_cube(int c_ndx, float *v) 00308 { 00309 m_case = cases[c_ndx][0]; 00310 m_config = cases[c_ndx][1]; 00311 m_subconfig = 0; 00312 00313 switch (m_case) { 00314 case 0: 00315 return -1; 00316 00317 case 1: 00318 return OFFSET_T1 + m_config; 00319 00320 case 2: 00321 return OFFSET_T2 + m_config; 00322 00323 case 3: 00324 if (mc33_test_face(test[OFFSET_TEST3 + m_config][0], v)) 00325 return OFFSET_T3_2 + m_config; /* 3.2 */ 00326 else 00327 return OFFSET_T3_1 + m_config; /* 3.1 */ 00328 00329 case 4: 00330 if (mc33_test_interior(test[OFFSET_TEST4 + m_config][0], v)) 00331 return OFFSET_T4_1 + m_config; /* 4.1 */ 00332 else 00333 return OFFSET_T4_2 + m_config; /* 4.2 */ 00334 00335 case 5: 00336 return OFFSET_T5 + m_config; 00337 00338 case 6: 00339 if (mc33_test_face(test[OFFSET_TEST6 + m_config][0], v)) 00340 return OFFSET_T6_2 + m_config; /* 6.2 */ 00341 else { 00342 if (mc33_test_interior(test[OFFSET_TEST6 + m_config][1], v)) 00343 return OFFSET_T6_1_1 + m_config; /* 6.1.1 */ 00344 else 00345 return OFFSET_T6_1_2 + m_config; /* 6.1.2 */ 00346 } 00347 00348 case 7: 00349 if (mc33_test_face(test[OFFSET_TEST7 + m_config][0], v)) 00350 m_subconfig += 1; 00351 if (mc33_test_face(test[OFFSET_TEST7 + m_config][1], v)) 00352 m_subconfig += 2; 00353 if (mc33_test_face(test[OFFSET_TEST7 + m_config][2], v)) 00354 m_subconfig += 4; 00355 00356 switch (subconfig7[m_subconfig]) { 00357 case 0: 00358 if (mc33_test_interior(test[OFFSET_TEST7 + m_config][3], v)) 00359 return OFFSET_T7_4_2 + m_config; /* 7.4.2 */ 00360 else 00361 return OFFSET_T7_4_1 + m_config; /* 7.4.1 */ 00362 case 1: 00363 return OFFSET_T7_3_S1 + m_config; /* 7.3 */ 00364 case 2: 00365 return OFFSET_T7_3_S2 + m_config; /* 7.3 */ 00366 case 3: 00367 return OFFSET_T7_3_S3 + m_config; /* 7.3 */ 00368 case 4: 00369 return OFFSET_T7_2_S1 + m_config; /* 7.2 */ 00370 case 5: 00371 return OFFSET_T7_2_S2 + m_config; /* 7.2 */ 00372 case 6: 00373 return OFFSET_T7_2_S3 + m_config; /* 7.2 */ 00374 case 7: 00375 return OFFSET_T7_1 + m_config; /* 7.1 */ 00376 }; 00377 00378 case 8: 00379 return OFFSET_T8 + m_config; 00380 00381 case 9: 00382 return OFFSET_T9 + m_config; 00383 00384 case 10: 00385 if (mc33_test_face(test[OFFSET_TEST10 + m_config][0], v)) { 00386 if (mc33_test_face(test[OFFSET_TEST10 + m_config][1], v)) 00387 return OFFSET_T10_1_1_S2 + m_config; /* 10.1.1 */ 00388 else { 00389 return OFFSET_T10_2_S1 + m_config; /* 10.2 */ 00390 } 00391 } 00392 else { 00393 if (mc33_test_face(test[OFFSET_TEST10 + m_config][1], v)) { 00394 return OFFSET_T10_2_S2 + m_config; /* 10.2 */ 00395 } 00396 else { 00397 if (mc33_test_interior(test[OFFSET_TEST10 + m_config][2], v)) 00398 return OFFSET_T10_1_1_S1 + m_config; /* 10.1.1 */ 00399 else 00400 return OFFSET_T10_1_2 + m_config; /* 10.1.2 */ 00401 } 00402 } 00403 00404 case 11: 00405 return OFFSET_T11 + m_config; 00406 00407 case 12: 00408 if (mc33_test_face(test[OFFSET_TEST12 + m_config][0], v)) { 00409 if (mc33_test_face(test[OFFSET_TEST12 + m_config][1], v)) 00410 return OFFSET_T12_1_1_S2 + m_config; /* 12.1.1 */ 00411 else { 00412 return OFFSET_T12_2_S1 + m_config; /* 12.2 */ 00413 } 00414 } 00415 else { 00416 if (mc33_test_face(test[OFFSET_TEST12 + m_config][1], v)) { 00417 return OFFSET_T12_2_S2 + m_config; /* 12.2 */ 00418 } 00419 else { 00420 if (mc33_test_interior(test[OFFSET_TEST12 + m_config][2], v)) 00421 return OFFSET_T12_1_1_S1 + m_config; /* 12.1.1 */ 00422 else 00423 return OFFSET_T12_1_2 + m_config; /* 12.1.2 */ 00424 } 00425 } 00426 00427 case 13: 00428 if (mc33_test_face(test[OFFSET_TEST13 + m_config][0], v)) 00429 m_subconfig += 1; 00430 if (mc33_test_face(test[OFFSET_TEST13 + m_config][1], v)) 00431 m_subconfig += 2; 00432 if (mc33_test_face(test[OFFSET_TEST13 + m_config][2], v)) 00433 m_subconfig += 4; 00434 if (mc33_test_face(test[OFFSET_TEST13 + m_config][3], v)) 00435 m_subconfig += 8; 00436 if (mc33_test_face(test[OFFSET_TEST13 + m_config][4], v)) 00437 m_subconfig += 16; 00438 if (mc33_test_face(test[OFFSET_TEST13 + m_config][5], v)) 00439 m_subconfig += 32; 00440 00441 switch (subconfig13[m_subconfig]) { 00442 case 0: /* 13.1 */ 00443 return OFFSET_T13_1_S1 + m_config; 00444 00445 case 1: /* 13.2 */ 00446 return OFFSET_T13_2_S1 + 0 + m_config * 6; 00447 case 2: /* 13.2 */ 00448 return OFFSET_T13_2_S1 + 1 + m_config * 6; 00449 case 3: /* 13.2 */ 00450 return OFFSET_T13_2_S1 + 2 + m_config * 6; 00451 case 4: /* 13.2 */ 00452 return OFFSET_T13_2_S1 + 3 + m_config * 6; 00453 case 5: /* 13.2 */ 00454 return OFFSET_T13_2_S1 + 4 + m_config * 6; 00455 case 6: /* 13.2 */ 00456 return OFFSET_T13_2_S1 + 5 + m_config * 6; 00457 00458 case 7: /* 13.3 */ 00459 return OFFSET_T13_3_S1 + 0 + m_config * 12; 00460 case 8: /* 13.3 */ 00461 return OFFSET_T13_3_S1 + 1 + m_config * 12; 00462 case 9: /* 13.3 */ 00463 return OFFSET_T13_3_S1 + 2 + m_config * 12; 00464 case 10: /* 13.3 */ 00465 return OFFSET_T13_3_S1 + 3 + m_config * 12; 00466 case 11: /* 13.3 */ 00467 return OFFSET_T13_3_S1 + 4 + m_config * 12; 00468 case 12: /* 13.3 */ 00469 return OFFSET_T13_3_S1 + 5 + m_config * 12; 00470 case 13: /* 13.3 */ 00471 return OFFSET_T13_3_S1 + 6 + m_config * 12; 00472 case 14: /* 13.3 */ 00473 return OFFSET_T13_3_S1 + 7 + m_config * 12; 00474 case 15: /* 13.3 */ 00475 return OFFSET_T13_3_S1 + 8 + m_config * 12; 00476 case 16: /* 13.3 */ 00477 return OFFSET_T13_3_S1 + 9 + m_config * 12; 00478 case 17: /* 13.3 */ 00479 return OFFSET_T13_3_S1 + 10 + m_config * 12; 00480 case 18: /* 13.3 */ 00481 return OFFSET_T13_3_S1 + 11 + m_config * 12; 00482 00483 case 19: /* 13.4 */ 00484 return OFFSET_T13_4 + 0 + m_config * 4; 00485 case 20: /* 13.4 */ 00486 return OFFSET_T13_4 + 1 + m_config * 4; 00487 case 21: /* 13.4 */ 00488 return OFFSET_T13_4 + 2 + m_config * 4; 00489 case 22: /* 13.4 */ 00490 return OFFSET_T13_4 + 3 + m_config * 4; 00491 00492 case 23: /* 13.5 */ 00493 m_subconfig = 0; 00494 if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v)) 00495 return OFFSET_T13_5_1 + 0 + m_config * 4; 00496 else 00497 return OFFSET_T13_5_2 + 0 + m_config * 4; 00498 case 24: /* 13.5 */ 00499 m_subconfig = 1; 00500 if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v)) 00501 return OFFSET_T13_5_1 + 1 + m_config * 4; 00502 else 00503 return OFFSET_T13_5_2 + 1 + m_config * 4; 00504 case 25: /* 13.5 */ 00505 m_subconfig = 2; 00506 if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v)) 00507 return OFFSET_T13_5_1 + 2 + m_config * 4; 00508 else 00509 return OFFSET_T13_5_2 + 2 + m_config * 4; 00510 case 26: /* 13.5 */ 00511 m_subconfig = 3; 00512 if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v)) 00513 return OFFSET_T13_5_1 + 3 + m_config * 4; 00514 else 00515 return OFFSET_T13_5_2 + 3 + m_config * 4; 00516 00517 case 27: /* 13.3 */ 00518 return OFFSET_T13_3_S2 + 0 + m_config * 12; 00519 case 28: /* 13.3 */ 00520 return OFFSET_T13_3_S2 + 1 + m_config * 12; 00521 case 29: /* 13.3 */ 00522 return OFFSET_T13_3_S2 + 2 + m_config * 12; 00523 case 30: /* 13.3 */ 00524 return OFFSET_T13_3_S2 + 3 + m_config * 12; 00525 case 31: /* 13.3 */ 00526 return OFFSET_T13_3_S2 + 4 + m_config * 12; 00527 case 32: /* 13.3 */ 00528 return OFFSET_T13_3_S2 + 5 + m_config * 12; 00529 case 33: /* 13.3 */ 00530 return OFFSET_T13_3_S2 + 6 + m_config * 12; 00531 case 34: /* 13.3 */ 00532 return OFFSET_T13_3_S2 + 7 + m_config * 12; 00533 case 35: /* 13.3 */ 00534 return OFFSET_T13_3_S2 + 8 + m_config * 12; 00535 case 36: /* 13.3 */ 00536 return OFFSET_T13_3_S2 + 9 + m_config * 12; 00537 case 37: /* 13.3 */ 00538 return OFFSET_T13_3_S2 + 10 + m_config * 12; 00539 case 38: /* 13.3 */ 00540 return OFFSET_T13_3_S2 + 11 + m_config * 12; 00541 00542 case 39: /* 13.2 */ 00543 return OFFSET_T13_2_S2 + 0 + m_config * 6; 00544 case 40: /* 13.2 */ 00545 return OFFSET_T13_2_S2 + 1 + m_config * 6; 00546 case 41: /* 13.2 */ 00547 return OFFSET_T13_2_S2 + 2 + m_config * 6; 00548 case 42: /* 13.2 */ 00549 return OFFSET_T13_2_S2 + 3 + m_config * 6; 00550 case 43: /* 13.2 */ 00551 return OFFSET_T13_2_S2 + 4 + m_config * 6; 00552 case 44: /* 13.2 */ 00553 return OFFSET_T13_2_S2 + 5 + m_config * 6; 00554 00555 case 45: /* 13.1 */ 00556 return OFFSET_T13_1_S2 + m_config; 00557 00558 default: 00559 fprintf(stderr, "Marching Cubes: Impossible case 13?\n"); 00560 } 00561 00562 case 14: 00563 return OFFSET_T14 + m_config; 00564 } 00565 00566 return -1; 00567 }