/* Based on implementation of MarchingCubes 33 Algorithm by Thomas Lewiner, thomas.lewiner@polytechnique.org, Math Dept, PUC-Rio Implementation in ogsf library Tomas Paudits, 2004 */ #include #include #include #include "mc33_table.h" unsigned char m_case, m_config, m_subconfig; /************************************************************************/ int mc33_test_face( char face , float *v) { float A,B,C,D ; switch( face ) { case -1 : case 1 : A = v[0] ; B = v[4] ; C = v[5] ; D = v[1] ; break ; case -2 : case 2 : A = v[1] ; B = v[5] ; C = v[6] ; D = v[2] ; break ; case -3 : case 3 : A = v[2] ; B = v[6] ; C = v[7] ; D = v[3] ; break ; case -4 : case 4 : A = v[3] ; B = v[7] ; C = v[4] ; D = v[0] ; break ; case -5 : case 5 : A = v[0] ; B = v[3] ; C = v[2] ; D = v[1] ; break ; case -6 : case 6 : A = v[4] ; B = v[7] ; C = v[6] ; D = v[5] ; break ; default : fprintf(stderr, "Invalid face code %d\n", face ); A = B = C = D = 0 ; }; return face * A * ( A*C - B*D ) >= 0; } /************************************************************************/ int mc33_test_interior( char s, float *v) { float t, At = 0, Bt = 0, Ct = 0, Dt = 0, a, b ; char test = 0 ; char edge = -1 ; switch( m_case ) { case 4 : case 10 : a = ( v[4] - v[0] ) * ( v[6] - v[2] ) - ( v[7] - v[3] ) * ( v[5] - v[1] ) ; b = v[2] * ( v[4] - v[0] ) + v[0] * ( v[6] - v[2] ) - v[1] * ( v[7] - v[3] ) - v[3] * ( v[5] - v[1] ) ; t = - b / (2*a) ; if( t<0 || t>1 ) return s > 0 ; At = v[0] + ( v[4] - v[0] ) * t ; Bt = v[3] + ( v[7] - v[3] ) * t ; Ct = v[2] + ( v[6] - v[2] ) * t ; Dt = v[1] + ( v[5] - v[1] ) * t ; break ; case 6 : case 7 : case 12 : case 13 : switch( m_case ) { case 6 : edge = cell_table[OFFSET_T6_1_1 + m_config].polys[0]; break ; case 7 : edge = cell_table[OFFSET_T7_4_1 + m_config].polys[13]; break ; case 12 : edge = cell_table[OFFSET_T12_2_S1 + m_config].polys[14]; break ; case 13 : edge = cell_table[OFFSET_T13_5_1 + m_subconfig + m_config * 4].polys[2]; break ; } switch( edge ) { case 0 : t = v[0] / ( v[0] - v[1] ) ; At = 0 ; Bt = v[3] + ( v[2] - v[3] ) * t ; Ct = v[7] + ( v[6] - v[7] ) * t ; Dt = v[4] + ( v[5] - v[4] ) * t ; break ; case 1 : t = v[1] / ( v[1] - v[2] ) ; At = 0 ; Bt = v[0] + ( v[3] - v[0] ) * t ; Ct = v[4] + ( v[7] - v[4] ) * t ; Dt = v[5] + ( v[6] - v[5] ) * t ; break ; case 2 : t = v[2] / ( v[2] - v[3] ) ; At = 0 ; Bt = v[1] + ( v[0] - v[1] ) * t ; Ct = v[5] + ( v[4] - v[5] ) * t ; Dt = v[6] + ( v[7] - v[6] ) * t ; break ; case 3 : t = v[3] / ( v[3] - v[0] ) ; At = 0 ; Bt = v[2] + ( v[1] - v[2] ) * t ; Ct = v[6] + ( v[5] - v[6] ) * t ; Dt = v[7] + ( v[4] - v[7] ) * t ; break ; case 4 : t = v[4] / ( v[4] - v[5] ) ; At = 0 ; Bt = v[7] + ( v[6] - v[7] ) * t ; Ct = v[3] + ( v[2] - v[3] ) * t ; Dt = v[0] + ( v[1] - v[0] ) * t ; break ; case 5 : t = v[5] / ( v[5] - v[6] ) ; At = 0 ; Bt = v[4] + ( v[7] - v[4] ) * t ; Ct = v[0] + ( v[3] - v[0] ) * t ; Dt = v[1] + ( v[2] - v[1] ) * t ; break ; case 6 : t = v[6] / ( v[6] - v[7] ) ; At = 0 ; Bt = v[5] + ( v[4] - v[5] ) * t ; Ct = v[1] + ( v[0] - v[1] ) * t ; Dt = v[2] + ( v[3] - v[2] ) * t ; break ; case 7 : t = v[7] / ( v[7] - v[4] ) ; At = 0 ; Bt = v[6] + ( v[5] - v[6] ) * t ; Ct = v[2] + ( v[1] - v[2] ) * t ; Dt = v[3] + ( v[0] - v[3] ) * t ; break ; case 8 : t = v[0] / ( v[0] - v[4] ) ; At = 0 ; Bt = v[3] + ( v[7] - v[3] ) * t ; Ct = v[2] + ( v[6] - v[2] ) * t ; Dt = v[1] + ( v[5] - v[1] ) * t ; break ; case 9 : t = v[1] / ( v[1] - v[5] ) ; At = 0 ; Bt = v[0] + ( v[4] - v[0] ) * t ; Ct = v[3] + ( v[7] - v[3] ) * t ; Dt = v[2] + ( v[6] - v[2] ) * t ; break ; case 10 : t = v[2] / ( v[2] - v[6] ) ; At = 0 ; Bt = v[1] + ( v[5] - v[1] ) * t ; Ct = v[0] + ( v[4] - v[0] ) * t ; Dt = v[3] + ( v[7] - v[3] ) * t ; break ; case 11 : t = v[3] / ( v[3] - v[7] ) ; At = 0 ; Bt = v[2] + ( v[6] - v[2] ) * t ; Ct = v[1] + ( v[5] - v[1] ) * t ; Dt = v[0] + ( v[4] - v[0] ) * t ; break ; default : fprintf(stderr, "Invalid edge %d\n", edge ) ; break ; } break ; default : fprintf(stderr, "Invalid ambiguous case %d\n", m_case ) ; break ; } if( At >= 0 ) test ++ ; if( Bt >= 0 ) test += 2 ; if( Ct >= 0 ) test += 4 ; if( Dt >= 0 ) test += 8 ; switch( test ) { case 0 : return s>0 ; case 1 : return s>0 ; case 2 : return s>0 ; case 3 : return s>0 ; case 4 : return s>0 ; case 5 : if( At * Ct < Bt * Dt ) return s>0 ; break ; case 6 : return s>0 ; case 7 : return s<0 ; case 8 : return s>0 ; case 9 : return s>0 ; case 10 : if( At * Ct >= Bt * Dt ) return s>0 ; break ; case 11 : return s<0 ; case 12 : return s>0 ; case 13 : return s<0 ; case 14 : return s<0 ; case 15 : return s<0 ; } return s < 0; } /************************************************************************/ int mc33_process_cube(int c_ndx, float *v) { m_case = cases[c_ndx][0] ; m_config = cases[c_ndx][1] ; m_subconfig = 0 ; switch( m_case ) { case 0 : return -1; case 1 : return OFFSET_T1 + m_config; case 2 : return OFFSET_T2 + m_config; case 3 : if (mc33_test_face(test[OFFSET_TEST3 + m_config][0], v)) return OFFSET_T3_2 + m_config; /* 3.2 */ else return OFFSET_T3_1 + m_config; /* 3.1 */ case 4 : if( mc33_test_interior(test[OFFSET_TEST4 + m_config][0], v)) return OFFSET_T4_1 + m_config; /* 4.1 */ else return OFFSET_T4_2 + m_config; /* 4.2 */ case 5 : return OFFSET_T5 + m_config; case 6 : if( mc33_test_face(test[OFFSET_TEST6 + m_config][0],v)) return OFFSET_T6_2 + m_config; /* 6.2 */ else { if( mc33_test_interior(test[OFFSET_TEST6 + m_config][1], v)) return OFFSET_T6_1_1 + m_config; /* 6.1.1 */ else return OFFSET_T6_1_2 + m_config; /* 6.1.2 */ } case 7 : if(mc33_test_face(test[OFFSET_TEST7 + m_config][0], v)) m_subconfig += 1 ; if(mc33_test_face(test[OFFSET_TEST7 + m_config][1], v)) m_subconfig += 2 ; if(mc33_test_face(test[OFFSET_TEST7 + m_config][2], v)) m_subconfig += 4 ; switch(subconfig7[m_subconfig]) { case 0 : if( mc33_test_interior(test[OFFSET_TEST7 + m_config][3], v)) return OFFSET_T7_4_2 + m_config; /* 7.4.2 */ else return OFFSET_T7_4_1 + m_config; /* 7.4.1 */ case 1 : return OFFSET_T7_3_S1 + m_config; /* 7.3 */ case 2 : return OFFSET_T7_3_S2 + m_config; /* 7.3 */ case 3 : return OFFSET_T7_3_S3 + m_config; /* 7.3 */ case 4 : return OFFSET_T7_2_S1 + m_config; /* 7.2 */ case 5 : return OFFSET_T7_2_S2 + m_config; /* 7.2 */ case 6 : return OFFSET_T7_2_S3 + m_config; /* 7.2 */ case 7 : return OFFSET_T7_1 + m_config; /* 7.1 */ }; case 8 : return OFFSET_T8 + m_config; case 9 : return OFFSET_T9 + m_config; case 10 : if(mc33_test_face(test[OFFSET_TEST10 + m_config][0], v)) { if(mc33_test_face(test[OFFSET_TEST10 + m_config][1], v)) return OFFSET_T10_1_1_S2 + m_config; /* 10.1.1 */ else { return OFFSET_T10_2_S1 + m_config; /* 10.2 */ } } else { if(mc33_test_face(test[OFFSET_TEST10 + m_config][1], v)) { return OFFSET_T10_2_S2 + m_config; /* 10.2 */ } else { if( mc33_test_interior(test[OFFSET_TEST10 + m_config][2], v)) return OFFSET_T10_1_1_S1 + m_config; /* 10.1.1 */ else return OFFSET_T10_1_2 + m_config; /* 10.1.2 */ } } case 11 : return OFFSET_T11 + m_config; case 12 : if(mc33_test_face(test[OFFSET_TEST12 + m_config][0], v)) { if(mc33_test_face(test[OFFSET_TEST12 + m_config][1], v)) return OFFSET_T12_1_1_S2 + m_config; /* 12.1.1 */ else { return OFFSET_T12_2_S1 + m_config; /* 12.2 */ } } else { if(mc33_test_face(test[OFFSET_TEST12 + m_config][1], v)) { return OFFSET_T12_2_S2 + m_config; /* 12.2 */ } else { if( mc33_test_interior(test[OFFSET_TEST12 + m_config][2], v)) return OFFSET_T12_1_1_S1 + m_config; /* 12.1.1 */ else return OFFSET_T12_1_2 + m_config; /* 12.1.2 */ } } case 13 : if(mc33_test_face(test[OFFSET_TEST13 + m_config][0], v)) m_subconfig += 1 ; if(mc33_test_face(test[OFFSET_TEST13 + m_config][1], v)) m_subconfig += 2 ; if(mc33_test_face(test[OFFSET_TEST13 + m_config][2], v)) m_subconfig += 4 ; if(mc33_test_face(test[OFFSET_TEST13 + m_config][3], v)) m_subconfig += 8 ; if(mc33_test_face(test[OFFSET_TEST13 + m_config][4], v)) m_subconfig += 16 ; if(mc33_test_face(test[OFFSET_TEST13 + m_config][5], v)) m_subconfig += 32 ; switch( subconfig13[m_subconfig] ) { case 0 :/* 13.1 */ return OFFSET_T13_1_S1 + m_config; case 1 :/* 13.2 */ return OFFSET_T13_2_S1 + 0 + m_config * 6; case 2 :/* 13.2 */ return OFFSET_T13_2_S1 + 1 + m_config * 6; case 3 :/* 13.2 */ return OFFSET_T13_2_S1 + 2 + m_config * 6; case 4 :/* 13.2 */ return OFFSET_T13_2_S1 + 3 + m_config * 6; case 5 :/* 13.2 */ return OFFSET_T13_2_S1 + 4 + m_config * 6; case 6 :/* 13.2 */ return OFFSET_T13_2_S1 + 5 + m_config * 6; case 7 :/* 13.3 */ return OFFSET_T13_3_S1 + 0 + m_config * 12; case 8 :/* 13.3 */ return OFFSET_T13_3_S1 + 1 + m_config * 12; case 9 :/* 13.3 */ return OFFSET_T13_3_S1 + 2 + m_config * 12; case 10 :/* 13.3 */ return OFFSET_T13_3_S1 + 3 + m_config * 12; case 11 :/* 13.3 */ return OFFSET_T13_3_S1 + 4 + m_config * 12; case 12 :/* 13.3 */ return OFFSET_T13_3_S1 + 5 + m_config * 12; case 13 :/* 13.3 */ return OFFSET_T13_3_S1 + 6 + m_config * 12; case 14 :/* 13.3 */ return OFFSET_T13_3_S1 + 7 + m_config * 12; case 15 :/* 13.3 */ return OFFSET_T13_3_S1 + 8 + m_config * 12; case 16 :/* 13.3 */ return OFFSET_T13_3_S1 + 9 + m_config * 12; case 17 :/* 13.3 */ return OFFSET_T13_3_S1 + 10 + m_config * 12; case 18 :/* 13.3 */ return OFFSET_T13_3_S1 + 11 + m_config * 12; case 19 :/* 13.4 */ return OFFSET_T13_4 + 0 + m_config * 4; case 20 :/* 13.4 */ return OFFSET_T13_4 + 1 + m_config * 4; case 21 :/* 13.4 */ return OFFSET_T13_4 + 2 + m_config * 4; case 22 :/* 13.4 */ return OFFSET_T13_4 + 3 + m_config * 4; case 23 :/* 13.5 */ m_subconfig = 0 ; if(mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v)) return OFFSET_T13_5_1 + 0 + m_config * 4; else return OFFSET_T13_5_2 + 0 + m_config * 4; case 24 :/* 13.5 */ m_subconfig = 1 ; if(mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v)) return OFFSET_T13_5_1 + 1 + m_config * 4; else return OFFSET_T13_5_2 + 1 + m_config * 4; case 25 :/* 13.5 */ m_subconfig = 2 ; if(mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v)) return OFFSET_T13_5_1 + 2 + m_config * 4; else return OFFSET_T13_5_2 + 2 + m_config * 4; case 26 :/* 13.5 */ m_subconfig = 3 ; if(mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v)) return OFFSET_T13_5_1 + 3 + m_config * 4; else return OFFSET_T13_5_2 + 3 + m_config * 4; case 27 :/* 13.3 */ return OFFSET_T13_3_S2 + 0 + m_config * 12; case 28 :/* 13.3 */ return OFFSET_T13_3_S2 + 1 + m_config * 12; case 29 :/* 13.3 */ return OFFSET_T13_3_S2 + 2 + m_config * 12; case 30 :/* 13.3 */ return OFFSET_T13_3_S2 + 3 + m_config * 12; case 31 :/* 13.3 */ return OFFSET_T13_3_S2 + 4 + m_config * 12; case 32 :/* 13.3 */ return OFFSET_T13_3_S2 + 5 + m_config * 12; case 33 :/* 13.3 */ return OFFSET_T13_3_S2 + 6 + m_config * 12; case 34 :/* 13.3 */ return OFFSET_T13_3_S2 + 7 + m_config * 12; case 35 :/* 13.3 */ return OFFSET_T13_3_S2 + 8 + m_config * 12; case 36 :/* 13.3 */ return OFFSET_T13_3_S2 + 9 + m_config * 12; case 37 :/* 13.3 */ return OFFSET_T13_3_S2 + 10 + m_config * 12; case 38 :/* 13.3 */ return OFFSET_T13_3_S2 + 11 + m_config * 12; case 39 :/* 13.2 */ return OFFSET_T13_2_S2 + 0 + m_config * 6; case 40 :/* 13.2 */ return OFFSET_T13_2_S2 + 1 + m_config * 6; case 41 :/* 13.2 */ return OFFSET_T13_2_S2 + 2 + m_config * 6; case 42 :/* 13.2 */ return OFFSET_T13_2_S2 + 3 + m_config * 6; case 43 :/* 13.2 */ return OFFSET_T13_2_S2 + 4 + m_config * 6; case 44 :/* 13.2 */ return OFFSET_T13_2_S2 + 5 + m_config * 6; case 45 :/* 13.1 */ return OFFSET_T13_1_S2 + m_config; default : fprintf(stderr, "Marching Cubes: Impossible case 13?\n" ) ; } case 14 : return OFFSET_T14 + m_config; } return -1; }