IsoExtractor.h (34407B)
1 // File : IsoExtractor.h 2 // Date : Fri Nov 25 19:12:52 2016 3 // Author : Fabian Wermelinger 4 // Description: Iso surface extractor 5 // Copyright 2016 ETH Zurich. All Rights Reserved. 6 #ifndef ISOEXTRACTOR_H_A1WQHFFC 7 #define ISOEXTRACTOR_H_A1WQHFFC 8 9 #include <vector> 10 #include <fstream> 11 #include <string> 12 #include <ctime> 13 #include <cstdlib> 14 #include <chrono> 15 #include "common.h" 16 #include "ArgumentParser.h" 17 #include "Interpolator.h" 18 19 struct Triangle 20 { 21 float vertices[3][3]; 22 } __attribute__((packed)); 23 24 static const Real gVertexOffset[8][3] = 25 { 26 {0.0, 0.0, 0.0},{1.0, 0.0, 0.0},{1.0, 1.0, 0.0},{0.0, 1.0, 0.0}, 27 {0.0, 0.0, 1.0},{1.0, 0.0, 1.0},{1.0, 1.0, 1.0},{0.0, 1.0, 1.0} 28 }; 29 30 static const int gEdgeConnection[12][2] = 31 { 32 {0,1}, {1,2}, {2,3}, {3,0}, 33 {4,5}, {5,6}, {6,7}, {7,4}, 34 {0,4}, {1,5}, {2,6}, {3,7} 35 }; 36 37 static const Real gEdgeDirection[12][3] = 38 { 39 {1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{-1.0, 0.0, 0.0},{0.0, -1.0, 0.0}, 40 {1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{-1.0, 0.0, 0.0},{0.0, -1.0, 0.0}, 41 {0.0, 0.0, 1.0},{0.0, 0.0, 1.0},{ 0.0, 0.0, 1.0},{0.0, 0.0, 1.0} 42 }; 43 44 static const int gTetrahedronEdgeConnection[6][2] = 45 { 46 {0,1}, {1,2}, {2,0}, {0,3}, {1,3}, {2,3} 47 }; 48 49 static const int gTetrahedronsInACube[6][4] = 50 { 51 {0,5,1,6}, 52 {0,1,2,6}, 53 {0,2,3,6}, 54 {0,3,7,6}, 55 {0,7,4,6}, 56 {0,4,5,6}, 57 }; 58 59 static const int gCubeEdgeFlags[256] = 60 { 61 0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, 62 0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, 63 0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435, 0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 64 0x3a0, 0x2a9, 0x1a3, 0x0aa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, 65 0x460, 0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60, 66 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, 67 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c, 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, 68 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0x0cc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 69 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0x0cc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0, 70 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x055, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650, 71 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5, 0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 72 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460, 73 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0, 74 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x033, 0x339, 0x230, 75 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190, 76 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x000 77 }; 78 79 static const int gTriangleConnectionTable[256][16] = 80 { 81 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 82 {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 83 {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 84 {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 85 {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 86 {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 87 {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 88 {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1}, 89 {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 90 {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 91 {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 92 {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1}, 93 {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 94 {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1}, 95 {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1}, 96 {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 97 {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 98 {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 99 {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 100 {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1}, 101 {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 102 {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1}, 103 {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, 104 {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1}, 105 {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 106 {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1}, 107 {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, 108 {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1}, 109 {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1}, 110 {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1}, 111 {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1}, 112 {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1}, 113 {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 114 {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 115 {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 116 {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1}, 117 {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 118 {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, 119 {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1}, 120 {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1}, 121 {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 122 {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, 123 {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, 124 {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1}, 125 {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1}, 126 {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1}, 127 {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1}, 128 {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1}, 129 {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 130 {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1}, 131 {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1}, 132 {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 133 {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1}, 134 {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1}, 135 {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1}, 136 {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1}, 137 {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1}, 138 {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1}, 139 {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1}, 140 {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1}, 141 {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1}, 142 {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1}, 143 {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1}, 144 {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 145 {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 146 {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 147 {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 148 {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, 149 {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 150 {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1}, 151 {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1}, 152 {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1}, 153 {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 154 {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, 155 {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, 156 {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1}, 157 {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1}, 158 {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1}, 159 {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1}, 160 {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1}, 161 {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 162 {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1}, 163 {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, 164 {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1}, 165 {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1}, 166 {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1}, 167 {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1}, 168 {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1}, 169 {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, 170 {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1}, 171 {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1}, 172 {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1}, 173 {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1}, 174 {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1}, 175 {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1}, 176 {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1}, 177 {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 178 {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1}, 179 {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1}, 180 {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1}, 181 {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1}, 182 {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1}, 183 {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 184 {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1}, 185 {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1}, 186 {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1}, 187 {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1}, 188 {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1}, 189 {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1}, 190 {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1}, 191 {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1}, 192 {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 193 {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1}, 194 {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1}, 195 {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1}, 196 {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1}, 197 {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1}, 198 {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1}, 199 {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1}, 200 {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 201 {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1}, 202 {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1}, 203 {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1}, 204 {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1}, 205 {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1}, 206 {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 207 {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1}, 208 {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 209 {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 210 {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 211 {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 212 {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, 213 {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 214 {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, 215 {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, 216 {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1}, 217 {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 218 {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1}, 219 {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1}, 220 {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1}, 221 {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1}, 222 {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1}, 223 {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1}, 224 {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1}, 225 {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 226 {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1}, 227 {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1}, 228 {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1}, 229 {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1}, 230 {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1}, 231 {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1}, 232 {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1}, 233 {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1}, 234 {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 235 {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1}, 236 {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1}, 237 {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1}, 238 {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1}, 239 {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1}, 240 {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 241 {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 242 {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, 243 {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, 244 {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1}, 245 {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, 246 {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1}, 247 {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1}, 248 {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1}, 249 {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1}, 250 {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1}, 251 {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1}, 252 {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1}, 253 {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1}, 254 {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1}, 255 {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1}, 256 {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1}, 257 {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1}, 258 {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1}, 259 {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1}, 260 {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1}, 261 {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1}, 262 {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1}, 263 {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1}, 264 {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1}, 265 {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1}, 266 {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1}, 267 {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1}, 268 {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 269 {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1}, 270 {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1}, 271 {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 272 {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 273 {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 274 {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1}, 275 {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1}, 276 {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1}, 277 {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1}, 278 {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1}, 279 {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1}, 280 {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1}, 281 {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1}, 282 {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1}, 283 {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1}, 284 {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1}, 285 {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 286 {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1}, 287 {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1}, 288 {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 289 {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1}, 290 {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1}, 291 {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1}, 292 {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1}, 293 {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1}, 294 {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1}, 295 {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1}, 296 {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 297 {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1}, 298 {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1}, 299 {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1}, 300 {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1}, 301 {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1}, 302 {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 303 {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1}, 304 {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 305 {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1}, 306 {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1}, 307 {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1}, 308 {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1}, 309 {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1}, 310 {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1}, 311 {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1}, 312 {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1}, 313 {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1}, 314 {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1}, 315 {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1}, 316 {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 317 {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1}, 318 {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1}, 319 {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 320 {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 321 {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 322 {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1}, 323 {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1}, 324 {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 325 {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1}, 326 {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1}, 327 {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 328 {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 329 {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1}, 330 {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 331 {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1}, 332 {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 333 {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 334 {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 335 {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 336 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1} 337 }; 338 339 static const int gTetrahedronEdgeFlags[16] = 340 { 341 0x00, 0x0d, 0x13, 0x1e, 0x26, 0x2b, 0x35, 0x38, 0x38, 0x35, 0x2b, 0x26, 0x1e, 0x13, 0x0d, 0x00 342 }; 343 344 static const int gTetrahedronTriangles[16][7] = 345 { 346 {-1, -1, -1, -1, -1, -1, -1}, 347 { 0, 3, 2, -1, -1, -1, -1}, 348 { 0, 1, 4, -1, -1, -1, -1}, 349 { 1, 4, 2, 2, 4, 3, -1}, 350 351 { 1, 2, 5, -1, -1, -1, -1}, 352 { 0, 3, 5, 0, 5, 1, -1}, 353 { 0, 2, 5, 0, 5, 4, -1}, 354 { 5, 4, 3, -1, -1, -1, -1}, 355 356 { 3, 4, 5, -1, -1, -1, -1}, 357 { 4, 5, 0, 5, 2, 0, -1}, 358 { 1, 5, 0, 5, 3, 0, -1}, 359 { 5, 2, 1, -1, -1, -1, -1}, 360 361 { 3, 4, 2, 2, 4, 1, -1}, 362 { 4, 1, 0, -1, -1, -1, -1}, 363 { 2, 3, 0, -1, -1, -1, -1}, 364 {-1, -1, -1, -1, -1, -1, -1} 365 }; 366 367 368 template <typename Tkernel, template<typename> class Tinterpolator=Interpolator> 369 class IsoExtractor 370 { 371 public: 372 IsoExtractor(ArgumentParser& parser) : 373 m_parser(parser), 374 m_interp(parser), 375 m_x0(parser("xmin").asDouble(0.0)), 376 m_y0(parser("ymin").asDouble(0.0)), 377 m_z0(parser("zmin").asDouble(0.0)) {} 378 ~IsoExtractor() {} 379 380 void extract(const Real isoval, const std::string fout) 381 { 382 const Real cscale = m_parser("cube_scale").asDouble(1.0); 383 assert(cscale > 0.0); 384 const int Nx = static_cast<Real>(m_interp.getNx()) / cscale; 385 const int Ny = static_cast<Real>(m_interp.getNy()) / cscale; 386 const int Nz = static_cast<Real>(m_interp.getNz()) / cscale; 387 assert(Nx > 0); 388 assert(Ny > 0); 389 assert(Nz > 0); 390 const Real chx = m_interp.getExtentX() / Nx; 391 const Real chy = m_interp.getExtentY() / Ny; 392 const Real chz = m_interp.getExtentZ() / Nz; 393 394 _extractIso(isoval, Nx, Ny, Nz, chx, chy, chz); 395 396 _writePLY(fout, m_triList); 397 } 398 399 400 protected: 401 ArgumentParser& m_parser; 402 Tinterpolator<Tkernel> m_interp; 403 const Real m_x0; 404 const Real m_y0; 405 const Real m_z0; 406 407 std::vector<Triangle> m_triList; 408 409 // helper 410 inline Real _getEdgeFrac(const Real v0, const Real v1, const Real t) 411 { 412 const Real dv = v1 - v0; 413 if (dv == 0.0) return 0.5; 414 return (t - v0)/dv; 415 } 416 417 inline void _extractIso(const Real isoval, 418 const int Nx, const int Ny, const int Nz, 419 const Real chx, const Real chy, const Real chz, 420 const Real Ox=0.0, const Real Oy=0.0, const Real Oz=0.0, 421 const bool verbose=true); 422 423 inline void _marchCube(const Real isoval, 424 const Real x, const Real y, const Real z, 425 const Real hx, const Real hy, const Real hz, 426 const Real Ox, const Real Oy, const Real Oz, 427 std::vector<Triangle>& tlist); 428 429 inline void _marchTetrahedron(const Real isoval, 430 const Real x, const Real y, const Real z, 431 const Real hx, const Real hy, const Real hz, 432 const Real Ox, const Real Oy, const Real Oz, 433 std::vector<Triangle>& tlist); 434 435 inline void _evalTetrahedron(const Real coords[4][3], const Real val[4], 436 const Real isoval, std::vector<Triangle>& tlist); 437 438 void _writePLY(const std::string basename, const std::vector<Triangle>& tList) 439 { 440 time_t rawtime; 441 std::time(&rawtime); 442 struct tm* timeinfo = std::localtime(&rawtime); 443 char buf[256]; 444 std::strftime(buf, 256, "%A, %h %d %Y, %r %Z %z", timeinfo); 445 446 const unsigned int one = 1; 447 const unsigned char* const lsb = (const unsigned char*)(&one); 448 const std::string endianess((*lsb) ? "binary_little_endian" : "binary_big_endian"); 449 450 // write header 451 std::ofstream ply(basename+".ply"); 452 ply << "ply" << std::endl; 453 ply << "format "; 454 if (m_parser("binary").asBool(true)) 455 ply << endianess; 456 else 457 ply << "ascii"; 458 ply << " 1.0" << std::endl; 459 ply << "comment generated by IsoExtractor.h on " << buf << std::endl; 460 ply << "element vertex " << tList.size()*3 << std::endl; 461 ply << "property float x" << std::endl; 462 ply << "property float y" << std::endl; 463 ply << "property float z" << std::endl; 464 ply << "element face " << tList.size() << std::endl; 465 ply << "property list uchar int vertex_index" << std::endl; 466 ply << "end_header" << std::endl; 467 ply.close(); 468 469 // write content 470 if (m_parser("binary").asBool(true)) 471 { 472 ply.open(basename+".ply", std::ios::out | std::ios::app | std::ios::binary); 473 _writePLY_binary(ply, tList); 474 } 475 else 476 { 477 ply.open(basename+".ply", std::ios::out | std::ios::app); 478 _writePLY_ascii(ply, tList); 479 } 480 ply.close(); 481 } 482 483 void _writePLY_binary(std::ostream& ply, const std::vector<Triangle>& tList) 484 { 485 for (auto t : tList) 486 { 487 float tc[3][3]; 488 for (int j = 0; j < 3; ++j) 489 { 490 tc[j][0] = t.vertices[j][0]; 491 tc[j][1] = t.vertices[j][1]; 492 tc[j][2] = t.vertices[j][2]; 493 } 494 ply.write((char*)&tc[0][0], 9*sizeof(float)); 495 } 496 497 const char nVert = 3; 498 for (int i = 0; i < tList.size()*3; i += 3) 499 { 500 int tv[3] = {i, i+1, i+2}; 501 ply.write(&nVert, sizeof(char)); 502 ply.write((char*)&tv[0], 3*sizeof(int)); 503 } 504 } 505 506 void _writePLY_ascii(std::ostream& ply, const std::vector<Triangle>& tList) 507 { 508 for (auto t : tList) 509 for (int j = 0; j < 3; ++j) 510 ply << t.vertices[j][0] << " " << t.vertices[j][1] << " " << t.vertices[j][2] << std::endl; 511 for (int i = 0; i < tList.size()*3; i += 3) 512 ply << "3 " << i << " " << i+1 << " " << i+2 << std::endl; 513 } 514 }; 515 516 517 template <typename Tkernel, template<typename> class Tinterpolator> 518 void IsoExtractor<Tkernel,Tinterpolator>::_extractIso(const Real isoval, 519 const int Nx, const int Ny, const int Nz, 520 const Real chx, const Real chy, const Real chz, 521 const Real Ox, const Real Oy, const Real Oz, const bool verbose) 522 { 523 const int speculate = m_parser("speculate_triangle").asInt(150000); 524 const std::string algorithm = m_parser("march").asString("cube"); 525 if (!(algorithm == "cube" || algorithm == "tetrahedron")) 526 { 527 if (verbose) 528 std::cerr << "ERROR: IsoExtractor: Unknown algorithm " << algorithm << std::endl; 529 abort(); 530 } 531 532 std::chrono::time_point<std::chrono::system_clock> start, end; 533 start = std::chrono::system_clock::now(); 534 if (verbose) 535 std::cout << "Extracting isovalue " << isoval << std::endl; 536 537 const size_t nCubes = static_cast<size_t>(Nx) * Ny * Nz; 538 size_t totalSize = 0; 539 #pragma omp parallel 540 { 541 std::vector<Triangle> myTriList; 542 myTriList.reserve(speculate); 543 544 if (algorithm == "cube") 545 { 546 #pragma omp for 547 for (size_t i = 0; i < nCubes; ++i) 548 { 549 const size_t ix = i%Nx; 550 const size_t iy = (i/Nx)%Ny; 551 const size_t iz = (i/(Nx*Ny))%Nz; 552 const Real x = ix*chx; 553 const Real y = iy*chy; 554 const Real z = iz*chz; 555 _marchCube(isoval, x, y, z, chx, chy, chz, Ox, Oy, Oz, myTriList); 556 } 557 } 558 else if (algorithm == "tetrahedron") 559 { 560 #pragma omp for 561 for (size_t i = 0; i < nCubes; ++i) 562 { 563 const size_t ix = i%Nx; 564 const size_t iy = (i/Nx)%Ny; 565 const size_t iz = (i/(Nx*Ny))%Nz; 566 const Real x = ix*chx; 567 const Real y = iy*chy; 568 const Real z = iz*chz; 569 _marchTetrahedron(isoval, x, y, z, chx, chy, chz, Ox, Oy, Oz, myTriList); 570 } 571 } 572 #pragma omp critical 573 { 574 totalSize += myTriList.size(); 575 } 576 #pragma omp master 577 { 578 m_triList.clear(); 579 m_triList.reserve(totalSize); 580 } 581 #pragma omp barrier 582 #pragma omp critical 583 { 584 m_triList.insert(m_triList.end(), myTriList.begin(), myTriList.end()); 585 } 586 } 587 end = std::chrono::system_clock::now(); 588 std::chrono::duration<double> elapsed_seconds = end-start; 589 if (verbose) 590 std::cout << "Extraction done. Elapsed time = " << elapsed_seconds.count() << "s" << std::endl; 591 } 592 593 594 template <typename Tkernel, template<typename> class Tinterpolator> 595 void IsoExtractor<Tkernel,Tinterpolator>::_marchCube(const Real isoval, 596 const Real x, const Real y, const Real z, 597 const Real hx, const Real hy, const Real hz, 598 const Real Ox, const Real Oy, const Real Oz, 599 std::vector<Triangle>& tlist) 600 { 601 Real valCubeVerts[8]; 602 Real coordEdgeVert[12][3]; 603 604 // compute cube vertex values 605 for (int i = 0; i < 8; ++i) 606 valCubeVerts[i] = m_interp( 607 x + hx*gVertexOffset[i][0], 608 y + hy*gVertexOffset[i][1], 609 z + hz*gVertexOffset[i][2]); 610 611 // find if there are vertices inside and outside the target value 612 int flagVert = 0; 613 for (int i = 0; i < 8; ++i) 614 { 615 if (valCubeVerts[i] <= isoval) flagVert |= 1 << i; 616 } 617 618 // indentify the edge intersections 619 const int flagEdge = gCubeEdgeFlags[flagVert]; 620 621 // if there are no intersections at all, bye 622 if(flagEdge == 0) 623 return; 624 625 // compute coordinates of triangle vertices on intersected edges 626 for (int i = 0; i < 12; ++i) 627 { 628 if (flagEdge & (1 << i)) 629 { 630 // linear interpolation 631 const Real frac = _getEdgeFrac(valCubeVerts[gEdgeConnection[i][0]], valCubeVerts[gEdgeConnection[i][1]], isoval); 632 coordEdgeVert[i][0] = x + hx*(gVertexOffset[gEdgeConnection[i][0]][0] + frac*gEdgeDirection[i][0]) + m_x0 + Ox; 633 coordEdgeVert[i][1] = y + hy*(gVertexOffset[gEdgeConnection[i][0]][1] + frac*gEdgeDirection[i][1]) + m_y0 + Oy; 634 coordEdgeVert[i][2] = z + hz*(gVertexOffset[gEdgeConnection[i][0]][2] + frac*gEdgeDirection[i][2]) + m_z0 + Oz; 635 } 636 } 637 638 // push the triangles found into the list 639 for (int i = 0; i < 5; ++i) 640 { 641 if (gTriangleConnectionTable[flagVert][3*i] < 0) break; 642 643 Triangle tri; 644 for (int j = 0; j < 3; ++j) 645 { 646 const int k = gTriangleConnectionTable[flagVert][3*i + j]; 647 tri.vertices[j][0] = coordEdgeVert[k][0]; 648 tri.vertices[j][1] = coordEdgeVert[k][1]; 649 tri.vertices[j][2] = coordEdgeVert[k][2]; 650 } 651 tlist.push_back(tri); 652 } 653 } 654 655 656 template <typename Tkernel, template<typename> class Tinterpolator> 657 void IsoExtractor<Tkernel,Tinterpolator>::_marchTetrahedron(const Real isoval, 658 const Real x, const Real y, const Real z, 659 const Real hx, const Real hy, const Real hz, 660 const Real Ox, const Real Oy, const Real Oz, 661 std::vector<Triangle>& tlist) 662 { 663 Real valCubeVerts[8]; 664 Real coordCubeVerts[8][3]; 665 Real valTetraVerts[4]; 666 Real coordTetraVerts[4][3]; 667 668 // compute cube vertex coordinates 669 for (int i = 0; i < 8; ++i) 670 { 671 coordCubeVerts[i][0] = x + hx*gVertexOffset[i][0]; 672 coordCubeVerts[i][1] = y + hy*gVertexOffset[i][1]; 673 coordCubeVerts[i][2] = z + hz*gVertexOffset[i][2]; 674 } 675 676 // compute cube vertex values 677 for (int i = 0; i < 8; ++i) 678 { 679 valCubeVerts[i] = m_interp( 680 coordCubeVerts[i][0], 681 coordCubeVerts[i][1], 682 coordCubeVerts[i][2]); 683 coordCubeVerts[i][0] += m_x0 + Ox; 684 coordCubeVerts[i][1] += m_y0 + Oy; 685 coordCubeVerts[i][2] += m_z0 + Oz; 686 } 687 688 for (int i = 0; i < 6; ++i) 689 { 690 for (int j = 0; j < 4; ++j) 691 { 692 const int vertID = gTetrahedronsInACube[i][j]; 693 coordTetraVerts[j][0] = coordCubeVerts[vertID][0]; 694 coordTetraVerts[j][1] = coordCubeVerts[vertID][1]; 695 coordTetraVerts[j][2] = coordCubeVerts[vertID][2]; 696 valTetraVerts[j] = valCubeVerts[vertID]; 697 } 698 _evalTetrahedron(coordTetraVerts, valTetraVerts, isoval, tlist); 699 } 700 } 701 702 703 template <typename Tkernel, template<typename> class Tinterpolator> 704 void IsoExtractor<Tkernel,Tinterpolator>::_evalTetrahedron(const Real coords[4][3], const Real val[4], 705 const Real isoval, 706 std::vector<Triangle>& tlist) 707 { 708 Real coordEdgeVert[6][3]; 709 710 // find if there are vertices inside and outside the target value 711 int flagVert = 0; 712 for (int i = 0; i < 4; ++i) 713 { 714 if (val[i] <= isoval) flagVert |= 1 << i; 715 } 716 717 // indentify the edge intersections 718 const int flagEdge = gTetrahedronEdgeFlags[flagVert]; 719 720 // if there are no intersections at all, bye 721 if(flagEdge == 0) 722 return; 723 724 // compute coordinates of triangle vertices on intersected edges 725 for (int i = 0; i < 6; ++i) 726 { 727 if (flagEdge & (1 << i)) 728 { 729 const int vert0 = gTetrahedronEdgeConnection[i][0]; 730 const int vert1 = gTetrahedronEdgeConnection[i][1]; 731 732 // linear interpolation 733 const Real frac = _getEdgeFrac(val[vert0], val[vert1], isoval); 734 const Real ifrac= 1.0 - frac; 735 736 coordEdgeVert[i][0] = ifrac*coords[vert0][0] + frac*coords[vert1][0]; 737 coordEdgeVert[i][1] = ifrac*coords[vert0][1] + frac*coords[vert1][1]; 738 coordEdgeVert[i][2] = ifrac*coords[vert0][2] + frac*coords[vert1][2]; 739 } 740 } 741 742 // push the triangles found into the list 743 for (int i = 0; i < 2; ++i) 744 { 745 if (gTetrahedronTriangles[flagVert][3*i] < 0) break; 746 747 Triangle tri; 748 for (int j = 0; j < 3; ++j) 749 { 750 const int k = gTetrahedronTriangles[flagVert][3*i + j]; 751 tri.vertices[j][0] = coordEdgeVert[k][0]; 752 tri.vertices[j][1] = coordEdgeVert[k][1]; 753 tri.vertices[j][2] = coordEdgeVert[k][2]; 754 } 755 tlist.push_back(tri); 756 } 757 } 758 759 #endif /* ISOEXTRACTOR_H_A1WQHFFC */