GRASS Programmer's Manual  6.4.2(2012)
gvl_calc2.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines