File indexing completed on 2024-04-28 03:40:49
0001 // Aqsis 0002 // Copyright (C) 2006, Paul C. Gregory 0003 // 0004 // Contact: pgregory@aqsis.org 0005 // 0006 // This library is free software; you can redistribute it and/or 0007 // modify it under the terms of the GNU General Public 0008 // License as published by the Free Software Foundation; either 0009 // version 2 of the License, or (at your option) any later version. 0010 // 0011 // This library is distributed in the hope that it will be useful, 0012 // but WITHOUT ANY WARRANTY; without even the implied warranty of 0013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0014 // General Public License for more details. 0015 // 0016 // You should have received a copy of the GNU General Public 0017 // License along with this library; if not, write to the Free Software 0018 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA 0019 0020 /** \file 0021 \brief MarchingCubes Algorithm 0022 \author Thomas Lewiner <thomas.lewiner@polytechnique.org> 0023 \author Math Dept, PUC-Rio 0024 \version 0.2 0025 \date 12/08/2002 0026 */ 0027 0028 // Minor modifications for KDE-Edu/Analitza Library: Copyright (C) 2014 by Percy Camilo T. Aucahuasi <percy.camilo.ta@gmail.com> 0029 0030 0031 #include "marchingcubes.h" 0032 0033 #include <math.h> 0034 #include <memory.h> 0035 #include <float.h> 0036 #include <stdio.h> 0037 0038 0039 0040 0041 #include "lookuptable.h" 0042 0043 // step size of the arrays of vertices and triangles 0044 static const int ALLOC_SIZE = 65536; 0045 0046 0047 //_____________________________________________________________________________ 0048 // print cube for debug 0049 void MarchingCubes::print_cube() 0050 { 0051 //Aqsis::log() << warning << i_cube[0] << " " << i_cube[1] << " " << i_cube[2] << " " << i_cube[3] << " " << i_cube[4] << " " << i_cube[5] << " " << i_cube[6] << " " << i_cube[7]) << std::endl; 0052 } 0053 0054 //_____________________________________________________________________________ 0055 // Constructor 0056 MarchingCubes::MarchingCubes( ) : 0057 i_originalMC(false), 0058 i_data ((double*)nullptr), 0059 i_x_verts (( int *)nullptr), 0060 i_y_verts (( int *)nullptr), 0061 i_z_verts (( int *)nullptr), 0062 i_nverts (0), 0063 i_ntrigs (0), 0064 i_Nverts (0), 0065 i_Ntrigs (0), 0066 i_vertices (( Vertex *)nullptr), 0067 i_triangles ((Triangle*)nullptr) 0068 { 0069 } 0070 //_____________________________________________________________________________ 0071 0072 void MarchingCubes::setupSpace(const SpaceLimits &spaceLimits) 0073 { 0074 /// 0075 clean_all(); 0076 0077 /// 0078 int x,y,z; 0079 int a = 64;// para mas de 10 de radio de mundo 80 para a es un buen valor 0080 x = a; 0081 y=a; 0082 z=a; 0083 0084 i_size_x =(x); 0085 i_size_y =(y); 0086 i_size_z =(z); 0087 0088 xmin = spaceLimits.minX; 0089 ymin = spaceLimits.minY; 0090 zmin = spaceLimits.minZ; 0091 xmax = spaceLimits.maxX; 0092 ymax = spaceLimits.maxY; 0093 zmax = spaceLimits.maxZ; 0094 0095 // qDebug() << xmin << xmax << "|" << ymin << ymax << "|" << zmin << zmax; 0096 0097 hx = (xmax-xmin)/i_size_x; 0098 hy = (ymax-ymin)/i_size_y; 0099 hz = (zmax-zmin)/i_size_z; 0100 0101 0102 /// 0103 init_all(); 0104 0105 } 0106 0107 0108 //_____________________________________________________________________________ 0109 // Destructor 0110 MarchingCubes::~MarchingCubes() 0111 //----------------------------------------------------------------------------- 0112 { 0113 clean_all() ; 0114 } 0115 //_____________________________________________________________________________ 0116 0117 0118 0119 //_____________________________________________________________________________ 0120 // main algorithm 0121 void MarchingCubes::run() 0122 //----------------------------------------------------------------------------- 0123 { 0124 /* 0125 long int tick = clock(); 0126 TqDouble perclocks = 1.0/(double)CLOCKS_PER_SEC; 0127 */ 0128 //pareciara que existen mas evialuaciones, pero no es asi ... en realidad 0129 //teniendo este bucle ahorr evaluaciones pues eval de analitza es muy costoso 0130 for( i_k = 0 ; i_k < i_size_z ; i_k++ ) 0131 for( i_j = 0 ; i_j < i_size_y ; i_j++ ) 0132 for( i_i = 0 ; i_i < i_size_x ; i_i++ ) 0133 { 0134 set_data (evalScalarField(xmin+hx*i_i, ymin+hy*i_j, zmin+hz*i_k), i_i, i_j, i_k); 0135 } 0136 0137 compute_intersection_points( ) ; 0138 0139 0140 0141 for( i_k = 0 ; i_k < i_size_z-1 ; i_k++ ) 0142 for( i_j = 0 ; i_j < i_size_y-1 ; i_j++ ) 0143 for( i_i = 0 ; i_i < i_size_x-1 ; i_i++ ) 0144 { 0145 i_lut_entry = 0 ; 0146 for( int p = 0 ; p < 8 ; ++p ) 0147 { 0148 i_cube[p] = get_data( i_i+((p^(p>>1))&1), i_j+((p>>1)&1), i_k+((p>>2)&1) ) ; 0149 if( fabs( i_cube[p] ) < FLT_EPSILON ) 0150 i_cube[p] = FLT_EPSILON ; 0151 if( i_cube[p] > 0 ) 0152 i_lut_entry += 1 << p ; 0153 } 0154 /* 0155 if( ( i_cube[0] = get_data( i_i , i_j , i_k ) ) > 0 ) i_lut_entry += 1 ; 0156 if( ( i_cube[1] = get_data(i_i+1, i_j , i_k ) ) > 0 ) i_lut_entry += 2 ; 0157 if( ( i_cube[2] = get_data(i_i+1,i_j+1, i_k ) ) > 0 ) i_lut_entry += 4 ; 0158 if( ( i_cube[3] = get_data( i_i ,i_j+1, i_k ) ) > 0 ) i_lut_entry += 8 ; 0159 if( ( i_cube[4] = get_data( i_i , i_j ,i_k+1) ) > 0 ) i_lut_entry += 16 ; 0160 if( ( i_cube[5] = get_data(i_i+1, i_j ,i_k+1) ) > 0 ) i_lut_entry += 32 ; 0161 if( ( i_cube[6] = get_data(i_i+1,i_j+1,i_k+1) ) > 0 ) i_lut_entry += 64 ; 0162 if( ( i_cube[7] = get_data(i_i ,i_j+1,i_k+1) ) > 0 ) i_lut_entry += 128 ; 0163 */ 0164 process_cube( ) ; 0165 } 0166 0167 /* 0168 Aqsis::log() << info << "the cpu tooks " << perclocks * (TqDouble) (clock() - tick ) << " secs." << std::endl; 0169 0170 for( i_i = 0 ; i_i < 15 ; i_i++ ) 0171 { 0172 Aqsis::log() << info << i_N[i_i] << " cases " << i_i << std::endl; 0173 } 0174 */ 0175 } 0176 //_____________________________________________________________________________ 0177 0178 0179 0180 //_____________________________________________________________________________ 0181 // init temporary structures (must set sizes before call) 0182 void MarchingCubes::init_temps() 0183 //----------------------------------------------------------------------------- 0184 { 0185 long int howmany = i_size_x * i_size_y * i_size_z; 0186 0187 i_data = new double[ howmany ]; 0188 i_x_verts = new int [ howmany ]; 0189 i_y_verts = new int [ howmany ]; 0190 i_z_verts = new int [ howmany ]; 0191 0192 while (!i_x_verts || !i_y_verts || !i_z_verts) 0193 { 0194 clean_temps(); 0195 i_size_x /= 2; 0196 i_size_y /= 2; 0197 i_size_z /= 2; 0198 howmany = i_size_x * i_size_y * i_size_z; 0199 0200 i_data = new double[ howmany ]; 0201 i_x_verts = new int [ howmany ]; 0202 i_y_verts = new int [ howmany ]; 0203 i_z_verts = new int [ howmany ]; 0204 0205 } 0206 memset( i_x_verts, -1, howmany * sizeof( int ) ) ; 0207 memset( i_y_verts, -1, howmany * sizeof( int ) ) ; 0208 memset( i_z_verts, -1, howmany * sizeof( int ) ) ; 0209 0210 memset( i_N, 0, 15 * sizeof(int) ) ; 0211 } 0212 //_____________________________________________________________________________ 0213 0214 0215 0216 //_____________________________________________________________________________ 0217 // init all structures (must set sizes before call) 0218 void MarchingCubes::init_all () 0219 //----------------------------------------------------------------------------- 0220 { 0221 init_temps() ; 0222 0223 i_nverts = i_ntrigs = 0 ; 0224 i_Nverts = i_Ntrigs = ALLOC_SIZE ; 0225 i_vertices = new Vertex [i_Nverts] ; 0226 i_triangles = new Triangle[i_Ntrigs] ; 0227 } 0228 //_____________________________________________________________________________ 0229 0230 0231 0232 //_____________________________________________________________________________ 0233 // clean temporary structures 0234 void MarchingCubes::clean_temps() 0235 //----------------------------------------------------------------------------- 0236 { 0237 if (i_data) 0238 delete [] i_data; 0239 if (i_x_verts) 0240 delete [] i_x_verts; 0241 if (i_y_verts) 0242 delete [] i_y_verts; 0243 if (i_z_verts) 0244 delete [] i_z_verts; 0245 0246 i_data = (double*)nullptr ; 0247 i_x_verts = (int*)nullptr ; 0248 i_y_verts = (int*)nullptr ; 0249 i_z_verts = (int*)nullptr ; 0250 } 0251 //_____________________________________________________________________________ 0252 0253 0254 0255 //_____________________________________________________________________________ 0256 // clean all structures 0257 void MarchingCubes::clean_all() 0258 //----------------------------------------------------------------------------- 0259 { 0260 clean_temps() ; 0261 delete [] i_vertices ; 0262 delete [] i_triangles ; 0263 i_vertices = (Vertex *)nullptr ; 0264 i_triangles = (Triangle *)nullptr ; 0265 i_nverts = i_ntrigs = 0 ; 0266 i_Nverts = i_Ntrigs = 0 ; 0267 0268 i_size_x = i_size_y = i_size_z = -1 ; 0269 } 0270 //_____________________________________________________________________________ 0271 0272 0273 0274 //_____________________________________________________________________________ 0275 //_____________________________________________________________________________ 0276 0277 0278 //_____________________________________________________________________________ 0279 // Compute the intersection points 0280 void MarchingCubes::compute_intersection_points( ) 0281 //----------------------------------------------------------------------------- 0282 { 0283 for( i_k = 0 ; i_k < i_size_z ; i_k++ ) 0284 for( i_j = 0 ; i_j < i_size_y ; i_j++ ) 0285 for( i_i = 0 ; i_i < i_size_x ; i_i++ ) 0286 { 0287 i_cube[0] = get_data( i_i, i_j, i_k ) ; 0288 if( i_i < i_size_x - 1 ) 0289 i_cube[1] = get_data(i_i+1, i_j , i_k ) ; 0290 else 0291 i_cube[1] = i_cube[0] ; 0292 0293 if( i_j < i_size_y - 1 ) 0294 i_cube[3] = get_data( i_i ,i_j+1, i_k ) ; 0295 else 0296 i_cube[3] = i_cube[0] ; 0297 0298 if( i_k < i_size_z - 1 ) 0299 i_cube[4] = get_data( i_i , i_j ,i_k+1) ; 0300 else 0301 i_cube[4] = i_cube[0] ; 0302 0303 if( fabs( i_cube[0] ) < FLT_EPSILON ) 0304 i_cube[0] = FLT_EPSILON ; 0305 if( fabs( i_cube[1] ) < FLT_EPSILON ) 0306 i_cube[1] = FLT_EPSILON ; 0307 if( fabs( i_cube[3] ) < FLT_EPSILON ) 0308 i_cube[3] = FLT_EPSILON ; 0309 if( fabs( i_cube[4] ) < FLT_EPSILON ) 0310 i_cube[4] = FLT_EPSILON ; 0311 0312 if( i_cube[0] < 0 ) 0313 { 0314 if( i_cube[1] > 0 ) 0315 set_x_vert( add_x_vertex( ), i_i,i_j,i_k ) ; 0316 if( i_cube[3] > 0 ) 0317 set_y_vert( add_y_vertex( ), i_i,i_j,i_k ) ; 0318 if( i_cube[4] > 0 ) 0319 set_z_vert( add_z_vertex( ), i_i,i_j,i_k ) ; 0320 } 0321 else 0322 { 0323 if( i_cube[1] < 0 ) 0324 set_x_vert( add_x_vertex( ), i_i,i_j,i_k ) ; 0325 if( i_cube[3] < 0 ) 0326 set_y_vert( add_y_vertex( ), i_i,i_j,i_k ) ; 0327 if( i_cube[4] < 0 ) 0328 set_z_vert( add_z_vertex( ), i_i,i_j,i_k ) ; 0329 } 0330 } 0331 } 0332 //_____________________________________________________________________________ 0333 0334 0335 0336 0337 0338 //_____________________________________________________________________________ 0339 // Test a face 0340 // if face>0 return true if the face contains a part of the surface 0341 bool MarchingCubes::test_face( int face ) 0342 //----------------------------------------------------------------------------- 0343 { 0344 double A,B,C,D ; 0345 0346 switch( face ) 0347 { 0348 case -1 : 0349 case 1 : 0350 A = i_cube[0] ; 0351 B = i_cube[4] ; 0352 C = i_cube[5] ; 0353 D = i_cube[1] ; 0354 break ; 0355 case -2 : 0356 case 2 : 0357 A = i_cube[1] ; 0358 B = i_cube[5] ; 0359 C = i_cube[6] ; 0360 D = i_cube[2] ; 0361 break ; 0362 case -3 : 0363 case 3 : 0364 A = i_cube[2] ; 0365 B = i_cube[6] ; 0366 C = i_cube[7] ; 0367 D = i_cube[3] ; 0368 break ; 0369 case -4 : 0370 case 4 : 0371 A = i_cube[3] ; 0372 B = i_cube[7] ; 0373 C = i_cube[4] ; 0374 D = i_cube[0] ; 0375 break ; 0376 case -5 : 0377 case 5 : 0378 A = i_cube[0] ; 0379 B = i_cube[3] ; 0380 C = i_cube[2] ; 0381 D = i_cube[1] ; 0382 break ; 0383 case -6 : 0384 case 6 : 0385 A = i_cube[4] ; 0386 B = i_cube[7] ; 0387 C = i_cube[6] ; 0388 D = i_cube[5] ; 0389 break ; 0390 default : 0391 // Aqsis::log() << warning << "Invalid face code " << face << std::endl; 0392 print_cube() ; 0393 A = B = C = D = 0 ; 0394 }; 0395 0396 return face * A * ( A*C - B*D ) >= 0 ; // face and A invert signs 0397 } 0398 //_____________________________________________________________________________ 0399 0400 0401 0402 0403 0404 //_____________________________________________________________________________ 0405 // Test the interior of a cube 0406 // if s == 7, return true if the interior is empty 0407 // if s ==-7, return false if the interior is empty 0408 bool MarchingCubes::test_interior( int s ) 0409 //----------------------------------------------------------------------------- 0410 { 0411 double t, At=0, Bt=0, Ct=0, Dt=0, a, b ; 0412 char test = 0 ; 0413 char edge = -1 ; // reference edge of the triangulation 0414 0415 switch( i_case ) 0416 { 0417 case 4 : 0418 case 10 : 0419 a = ( i_cube[4] - i_cube[0] ) * ( i_cube[6] - i_cube[2] ) - ( i_cube[7] - i_cube[3] ) * ( i_cube[5] - i_cube[1] ) ; 0420 b = i_cube[2] * ( i_cube[4] - i_cube[0] ) + i_cube[0] * ( i_cube[6] - i_cube[2] ) 0421 - i_cube[1] * ( i_cube[7] - i_cube[3] ) - i_cube[3] * ( i_cube[5] - i_cube[1] ) ; 0422 t = - b / (2*a) ; 0423 if( t<0 || t>1 ) 0424 return s>0 ; 0425 0426 At = i_cube[0] + ( i_cube[4] - i_cube[0] ) * t ; 0427 Bt = i_cube[3] + ( i_cube[7] - i_cube[3] ) * t ; 0428 Ct = i_cube[2] + ( i_cube[6] - i_cube[2] ) * t ; 0429 Dt = i_cube[1] + ( i_cube[5] - i_cube[1] ) * t ; 0430 break ; 0431 0432 case 6 : 0433 case 7 : 0434 case 12 : 0435 case 13 : 0436 switch( i_case ) 0437 { 0438 case 6 : 0439 edge = test6 [i_config][2] ; 0440 break ; 0441 case 7 : 0442 edge = test7 [i_config][4] ; 0443 break ; 0444 case 12 : 0445 edge = test12[i_config][3] ; 0446 break ; 0447 case 13 : 0448 edge = tiling13_5_1[i_config][i_subconfig][0] ; 0449 break ; 0450 } 0451 switch( edge ) 0452 { 0453 case 0 : 0454 t = i_cube[0] / ( i_cube[0] - i_cube[1] ) ; 0455 At = 0 ; 0456 Bt = i_cube[3] + ( i_cube[2] - i_cube[3] ) * t ; 0457 Ct = i_cube[7] + ( i_cube[6] - i_cube[7] ) * t ; 0458 Dt = i_cube[4] + ( i_cube[5] - i_cube[4] ) * t ; 0459 break ; 0460 case 1 : 0461 t = i_cube[1] / ( i_cube[1] - i_cube[2] ) ; 0462 At = 0 ; 0463 Bt = i_cube[0] + ( i_cube[3] - i_cube[0] ) * t ; 0464 Ct = i_cube[4] + ( i_cube[7] - i_cube[4] ) * t ; 0465 Dt = i_cube[5] + ( i_cube[6] - i_cube[5] ) * t ; 0466 break ; 0467 case 2 : 0468 t = i_cube[2] / ( i_cube[2] - i_cube[3] ) ; 0469 At = 0 ; 0470 Bt = i_cube[1] + ( i_cube[0] - i_cube[1] ) * t ; 0471 Ct = i_cube[5] + ( i_cube[4] - i_cube[5] ) * t ; 0472 Dt = i_cube[6] + ( i_cube[7] - i_cube[6] ) * t ; 0473 break ; 0474 case 3 : 0475 t = i_cube[3] / ( i_cube[3] - i_cube[0] ) ; 0476 At = 0 ; 0477 Bt = i_cube[2] + ( i_cube[1] - i_cube[2] ) * t ; 0478 Ct = i_cube[6] + ( i_cube[5] - i_cube[6] ) * t ; 0479 Dt = i_cube[7] + ( i_cube[4] - i_cube[7] ) * t ; 0480 break ; 0481 case 4 : 0482 t = i_cube[4] / ( i_cube[4] - i_cube[5] ) ; 0483 At = 0 ; 0484 Bt = i_cube[7] + ( i_cube[6] - i_cube[7] ) * t ; 0485 Ct = i_cube[3] + ( i_cube[2] - i_cube[3] ) * t ; 0486 Dt = i_cube[0] + ( i_cube[1] - i_cube[0] ) * t ; 0487 break ; 0488 case 5 : 0489 t = i_cube[5] / ( i_cube[5] - i_cube[6] ) ; 0490 At = 0 ; 0491 Bt = i_cube[4] + ( i_cube[7] - i_cube[4] ) * t ; 0492 Ct = i_cube[0] + ( i_cube[3] - i_cube[0] ) * t ; 0493 Dt = i_cube[1] + ( i_cube[2] - i_cube[1] ) * t ; 0494 break ; 0495 case 6 : 0496 t = i_cube[6] / ( i_cube[6] - i_cube[7] ) ; 0497 At = 0 ; 0498 Bt = i_cube[5] + ( i_cube[4] - i_cube[5] ) * t ; 0499 Ct = i_cube[1] + ( i_cube[0] - i_cube[1] ) * t ; 0500 Dt = i_cube[2] + ( i_cube[3] - i_cube[2] ) * t ; 0501 break ; 0502 case 7 : 0503 t = i_cube[7] / ( i_cube[7] - i_cube[4] ) ; 0504 At = 0 ; 0505 Bt = i_cube[6] + ( i_cube[5] - i_cube[6] ) * t ; 0506 Ct = i_cube[2] + ( i_cube[1] - i_cube[2] ) * t ; 0507 Dt = i_cube[3] + ( i_cube[0] - i_cube[3] ) * t ; 0508 break ; 0509 case 8 : 0510 t = i_cube[0] / ( i_cube[0] - i_cube[4] ) ; 0511 At = 0 ; 0512 Bt = i_cube[3] + ( i_cube[7] - i_cube[3] ) * t ; 0513 Ct = i_cube[2] + ( i_cube[6] - i_cube[2] ) * t ; 0514 Dt = i_cube[1] + ( i_cube[5] - i_cube[1] ) * t ; 0515 break ; 0516 case 9 : 0517 t = i_cube[1] / ( i_cube[1] - i_cube[5] ) ; 0518 At = 0 ; 0519 Bt = i_cube[0] + ( i_cube[4] - i_cube[0] ) * t ; 0520 Ct = i_cube[3] + ( i_cube[7] - i_cube[3] ) * t ; 0521 Dt = i_cube[2] + ( i_cube[6] - i_cube[2] ) * t ; 0522 break ; 0523 case 10 : 0524 t = i_cube[2] / ( i_cube[2] - i_cube[6] ) ; 0525 At = 0 ; 0526 Bt = i_cube[1] + ( i_cube[5] - i_cube[1] ) * t ; 0527 Ct = i_cube[0] + ( i_cube[4] - i_cube[0] ) * t ; 0528 Dt = i_cube[3] + ( i_cube[7] - i_cube[3] ) * t ; 0529 break ; 0530 case 11 : 0531 t = i_cube[3] / ( i_cube[3] - i_cube[7] ) ; 0532 At = 0 ; 0533 Bt = i_cube[2] + ( i_cube[6] - i_cube[2] ) * t ; 0534 Ct = i_cube[1] + ( i_cube[5] - i_cube[1] ) * t ; 0535 Dt = i_cube[0] + ( i_cube[4] - i_cube[0] ) * t ; 0536 break ; 0537 default : 0538 // Aqsis::log() << warning << "Invalid edge " << edge << std::endl; 0539 print_cube() ; 0540 break ; 0541 } 0542 break ; 0543 0544 default : 0545 // Aqsis::log() << warning << "invalid ambiguous case " << i_case << std::endl; 0546 print_cube() ; 0547 break ; 0548 } 0549 0550 if( At >= 0 ) 0551 test ++ ; 0552 if( Bt >= 0 ) 0553 test += 2 ; 0554 if( Ct >= 0 ) 0555 test += 4 ; 0556 if( Dt >= 0 ) 0557 test += 8 ; 0558 switch( test ) 0559 { 0560 case 0 : 0561 return s>0 ; 0562 case 1 : 0563 return s>0 ; 0564 case 2 : 0565 return s>0 ; 0566 case 3 : 0567 return s>0 ; 0568 case 4 : 0569 return s>0 ; 0570 case 5 : 0571 if( At * Ct < Bt * Dt ) 0572 return s>0 ; 0573 break ; 0574 case 6 : 0575 return s>0 ; 0576 case 7 : 0577 return s<0 ; 0578 case 8 : 0579 return s>0 ; 0580 case 9 : 0581 return s>0 ; 0582 case 10 : 0583 if( At * Ct >= Bt * Dt ) 0584 return s>0 ; 0585 break ; 0586 case 11 : 0587 return s<0 ; 0588 case 12 : 0589 return s>0 ; 0590 case 13 : 0591 return s<0 ; 0592 case 14 : 0593 return s<0 ; 0594 case 15 : 0595 return s<0 ; 0596 } 0597 0598 return s<0 ; 0599 } 0600 //_____________________________________________________________________________ 0601 0602 0603 0604 0605 //_____________________________________________________________________________ 0606 // Process a unit cube 0607 void MarchingCubes::process_cube( ) 0608 //----------------------------------------------------------------------------- 0609 { 0610 if( i_originalMC ) 0611 { 0612 char nt = 0 ; 0613 while( casesClassic[i_lut_entry][3*nt] != -1 ) 0614 nt++ ; 0615 add_triangle( casesClassic[i_lut_entry], nt ) ; 0616 return ; 0617 } 0618 0619 int v12 = -1 ; 0620 i_case = cases[i_lut_entry][0] ; 0621 i_config = cases[i_lut_entry][1] ; 0622 i_subconfig = 0 ; 0623 0624 i_N[i_case]++ ; 0625 0626 switch( i_case ) 0627 { 0628 case 0 : 0629 break ; 0630 0631 case 1 : 0632 add_triangle( tiling1[i_config], 1 ) ; 0633 break ; 0634 0635 case 2 : 0636 add_triangle( tiling2[i_config], 2 ) ; 0637 break ; 0638 0639 case 3 : 0640 if( test_face( test3[i_config]) ) 0641 add_triangle( tiling3_2[i_config], 4 ) ; // 3.2 0642 else 0643 add_triangle( tiling3_1[i_config], 2 ) ; // 3.1 0644 break ; 0645 0646 case 4 : 0647 if( test_interior( test4[i_config]) ) 0648 add_triangle( tiling4_1[i_config], 2 ) ; // 4.1.1 0649 else 0650 add_triangle( tiling4_2[i_config], 6 ) ; // 4.1.2 0651 break ; 0652 0653 case 5 : 0654 add_triangle( tiling5[i_config], 3 ) ; 0655 break ; 0656 0657 case 6 : 0658 if( test_face( test6[i_config][0]) ) 0659 add_triangle( tiling6_2[i_config], 5 ) ; // 6.2 0660 else 0661 { 0662 if( test_interior( test6[i_config][1]) ) 0663 add_triangle( tiling6_1_1[i_config], 3 ) ; // 6.1.1 0664 else 0665 add_triangle( tiling6_1_2[i_config], 7 ) ; // 6.1.2 0666 } 0667 break ; 0668 0669 case 7 : 0670 if( test_face( test7[i_config][0] ) ) 0671 i_subconfig += 1 ; 0672 if( test_face( test7[i_config][1] ) ) 0673 i_subconfig += 2 ; 0674 if( test_face( test7[i_config][2] ) ) 0675 i_subconfig += 4 ; 0676 switch( i_subconfig ) 0677 { 0678 case 0 : 0679 add_triangle( tiling7_1[i_config], 3 ) ; 0680 break ; 0681 case 1 : 0682 add_triangle( tiling7_2[i_config][0], 5 ) ; 0683 break ; 0684 case 2 : 0685 add_triangle( tiling7_2[i_config][1], 5 ) ; 0686 break ; 0687 case 3 : 0688 v12 = add_c_vertex() ; 0689 add_triangle( tiling7_3[i_config][0], 9, v12 ) ; 0690 break ; 0691 case 4 : 0692 add_triangle( tiling7_2[i_config][2], 5 ) ; 0693 break ; 0694 case 5 : 0695 v12 = add_c_vertex() ; 0696 add_triangle( tiling7_3[i_config][1], 9, v12 ) ; 0697 break ; 0698 case 6 : 0699 v12 = add_c_vertex() ; 0700 add_triangle( tiling7_3[i_config][2], 9, v12 ) ; 0701 break ; 0702 case 7 : 0703 if( test_interior( test7[i_config][3]) ) 0704 add_triangle( tiling7_4_2[i_config], 9 ) ; 0705 else 0706 add_triangle( tiling7_4_1[i_config], 5 ) ; 0707 break ; 0708 }; 0709 break ; 0710 0711 case 8 : 0712 add_triangle( tiling8[i_config], 2 ) ; 0713 break ; 0714 0715 case 9 : 0716 add_triangle( tiling9[i_config], 4 ) ; 0717 break ; 0718 0719 case 10 : 0720 if( test_face( test10[i_config][0]) ) 0721 { 0722 if( test_face( test10[i_config][1]) ) 0723 add_triangle( tiling10_1_1_[i_config], 4 ) ; // 10.1.1 0724 else 0725 { 0726 v12 = add_c_vertex() ; 0727 add_triangle( tiling10_2[i_config], 8, v12 ) ; // 10.2 0728 } 0729 } 0730 else 0731 { 0732 if( test_face( test10[i_config][1]) ) 0733 { 0734 v12 = add_c_vertex() ; 0735 add_triangle( tiling10_2_[i_config], 8, v12 ) ; // 10.2 0736 } 0737 else 0738 { 0739 if( test_interior( test10[i_config][2]) ) 0740 add_triangle( tiling10_1_1[i_config], 4 ) ; // 10.1.1 0741 else 0742 add_triangle( tiling10_1_2[i_config], 8 ) ; // 10.1.2 0743 } 0744 } 0745 break ; 0746 0747 case 11 : 0748 add_triangle( tiling11[i_config], 4 ) ; 0749 break ; 0750 0751 case 12 : 0752 if( test_face( test12[i_config][0]) ) 0753 { 0754 if( test_face( test12[i_config][1]) ) 0755 add_triangle( tiling12_1_1_[i_config], 4 ) ; // 12.1.1 0756 else 0757 { 0758 v12 = add_c_vertex() ; 0759 add_triangle( tiling12_2[i_config], 8, v12 ) ; // 12.2 0760 } 0761 } 0762 else 0763 { 0764 if( test_face( test12[i_config][1]) ) 0765 { 0766 v12 = add_c_vertex() ; 0767 add_triangle( tiling12_2_[i_config], 8, v12 ) ; // 12.2 0768 } 0769 else 0770 { 0771 if( test_interior( test12[i_config][2]) ) 0772 add_triangle( tiling12_1_1[i_config], 4 ) ; // 12.1.1 0773 else 0774 add_triangle( tiling12_1_2[i_config], 8 ) ; // 12.1.2 0775 } 0776 } 0777 break ; 0778 0779 case 13 : 0780 if( test_face( test13[i_config][0] ) ) 0781 i_subconfig += 1 ; 0782 if( test_face( test13[i_config][1] ) ) 0783 i_subconfig += 2 ; 0784 if( test_face( test13[i_config][2] ) ) 0785 i_subconfig += 4 ; 0786 if( test_face( test13[i_config][3] ) ) 0787 i_subconfig += 8 ; 0788 if( test_face( test13[i_config][4] ) ) 0789 i_subconfig += 16 ; 0790 if( test_face( test13[i_config][5] ) ) 0791 i_subconfig += 32 ; 0792 switch( subconfig13[i_subconfig] ) 0793 { 0794 case 0 :/* 13.1 */ 0795 add_triangle( tiling13_1[i_config], 4 ) ; 0796 break ; 0797 0798 case 1 :/* 13.2 */ 0799 add_triangle( tiling13_2[i_config][0], 6 ) ; 0800 break ; 0801 case 2 :/* 13.2 */ 0802 add_triangle( tiling13_2[i_config][1], 6 ) ; 0803 break ; 0804 case 3 :/* 13.2 */ 0805 add_triangle( tiling13_2[i_config][2], 6 ) ; 0806 break ; 0807 case 4 :/* 13.2 */ 0808 add_triangle( tiling13_2[i_config][3], 6 ) ; 0809 break ; 0810 case 5 :/* 13.2 */ 0811 add_triangle( tiling13_2[i_config][4], 6 ) ; 0812 break ; 0813 case 6 :/* 13.2 */ 0814 add_triangle( tiling13_2[i_config][5], 6 ) ; 0815 break ; 0816 0817 case 7 :/* 13.3 */ 0818 v12 = add_c_vertex() ; 0819 add_triangle( tiling13_3[i_config][0], 10, v12 ) ; 0820 break ; 0821 case 8 :/* 13.3 */ 0822 v12 = add_c_vertex() ; 0823 add_triangle( tiling13_3[i_config][1], 10, v12 ) ; 0824 break ; 0825 case 9 :/* 13.3 */ 0826 v12 = add_c_vertex() ; 0827 add_triangle( tiling13_3[i_config][2], 10, v12 ) ; 0828 break ; 0829 case 10 :/* 13.3 */ 0830 v12 = add_c_vertex() ; 0831 add_triangle( tiling13_3[i_config][3], 10, v12 ) ; 0832 break ; 0833 case 11 :/* 13.3 */ 0834 v12 = add_c_vertex() ; 0835 add_triangle( tiling13_3[i_config][4], 10, v12 ) ; 0836 break ; 0837 case 12 :/* 13.3 */ 0838 v12 = add_c_vertex() ; 0839 add_triangle( tiling13_3[i_config][5], 10, v12 ) ; 0840 break ; 0841 case 13 :/* 13.3 */ 0842 v12 = add_c_vertex() ; 0843 add_triangle( tiling13_3[i_config][6], 10, v12 ) ; 0844 break ; 0845 case 14 :/* 13.3 */ 0846 v12 = add_c_vertex() ; 0847 add_triangle( tiling13_3[i_config][7], 10, v12 ) ; 0848 break ; 0849 case 15 :/* 13.3 */ 0850 v12 = add_c_vertex() ; 0851 add_triangle( tiling13_3[i_config][8], 10, v12 ) ; 0852 break ; 0853 case 16 :/* 13.3 */ 0854 v12 = add_c_vertex() ; 0855 add_triangle( tiling13_3[i_config][9], 10, v12 ) ; 0856 break ; 0857 case 17 :/* 13.3 */ 0858 v12 = add_c_vertex() ; 0859 add_triangle( tiling13_3[i_config][10], 10, v12 ) ; 0860 break ; 0861 case 18 :/* 13.3 */ 0862 v12 = add_c_vertex() ; 0863 add_triangle( tiling13_3[i_config][11], 10, v12 ) ; 0864 break ; 0865 0866 case 19 :/* 13.4 */ 0867 v12 = add_c_vertex() ; 0868 add_triangle( tiling13_4[i_config][0], 12, v12 ) ; 0869 break ; 0870 case 20 :/* 13.4 */ 0871 v12 = add_c_vertex() ; 0872 add_triangle( tiling13_4[i_config][1], 12, v12 ) ; 0873 break ; 0874 case 21 :/* 13.4 */ 0875 v12 = add_c_vertex() ; 0876 add_triangle( tiling13_4[i_config][2], 12, v12 ) ; 0877 break ; 0878 case 22 :/* 13.4 */ 0879 v12 = add_c_vertex() ; 0880 add_triangle( tiling13_4[i_config][3], 12, v12 ) ; 0881 break ; 0882 0883 case 23 :/* 13.5 */ 0884 i_subconfig = 0 ; 0885 if( test_interior( test13[i_config][6] ) ) 0886 add_triangle( tiling13_5_1[i_config][0], 6 ) ; 0887 else 0888 add_triangle( tiling13_5_2[i_config][0], 10 ) ; 0889 break ; 0890 case 24 :/* 13.5 */ 0891 i_subconfig = 1 ; 0892 if( test_interior( test13[i_config][6] ) ) 0893 add_triangle( tiling13_5_1[i_config][1], 6 ) ; 0894 else 0895 add_triangle( tiling13_5_2[i_config][1], 10 ) ; 0896 break ; 0897 case 25 :/* 13.5 */ 0898 i_subconfig = 2 ; 0899 if( test_interior( test13[i_config][6] ) ) 0900 add_triangle( tiling13_5_1[i_config][2], 6 ) ; 0901 else 0902 add_triangle( tiling13_5_2[i_config][2], 10 ) ; 0903 break ; 0904 case 26 :/* 13.5 */ 0905 i_subconfig = 3 ; 0906 if( test_interior( test13[i_config][6] ) ) 0907 add_triangle( tiling13_5_1[i_config][3], 6 ) ; 0908 else 0909 add_triangle( tiling13_5_2[i_config][3], 10 ) ; 0910 break ; 0911 0912 case 27 :/* 13.3 */ 0913 v12 = add_c_vertex() ; 0914 add_triangle( tiling13_3_[i_config][0], 10, v12 ) ; 0915 break ; 0916 case 28 :/* 13.3 */ 0917 v12 = add_c_vertex() ; 0918 add_triangle( tiling13_3_[i_config][1], 10, v12 ) ; 0919 break ; 0920 case 29 :/* 13.3 */ 0921 v12 = add_c_vertex() ; 0922 add_triangle( tiling13_3_[i_config][2], 10, v12 ) ; 0923 break ; 0924 case 30 :/* 13.3 */ 0925 v12 = add_c_vertex() ; 0926 add_triangle( tiling13_3_[i_config][3], 10, v12 ) ; 0927 break ; 0928 case 31 :/* 13.3 */ 0929 v12 = add_c_vertex() ; 0930 add_triangle( tiling13_3_[i_config][4], 10, v12 ) ; 0931 break ; 0932 case 32 :/* 13.3 */ 0933 v12 = add_c_vertex() ; 0934 add_triangle( tiling13_3_[i_config][5], 10, v12 ) ; 0935 break ; 0936 case 33 :/* 13.3 */ 0937 v12 = add_c_vertex() ; 0938 add_triangle( tiling13_3_[i_config][6], 10, v12 ) ; 0939 break ; 0940 case 34 :/* 13.3 */ 0941 v12 = add_c_vertex() ; 0942 add_triangle( tiling13_3_[i_config][7], 10, v12 ) ; 0943 break ; 0944 case 35 :/* 13.3 */ 0945 v12 = add_c_vertex() ; 0946 add_triangle( tiling13_3_[i_config][8], 10, v12 ) ; 0947 break ; 0948 case 36 :/* 13.3 */ 0949 v12 = add_c_vertex() ; 0950 add_triangle( tiling13_3_[i_config][9], 10, v12 ) ; 0951 break ; 0952 case 37 :/* 13.3 */ 0953 v12 = add_c_vertex() ; 0954 add_triangle( tiling13_3_[i_config][10], 10, v12 ) ; 0955 break ; 0956 case 38 :/* 13.3 */ 0957 v12 = add_c_vertex() ; 0958 add_triangle( tiling13_3_[i_config][11], 10, v12 ) ; 0959 break ; 0960 0961 case 39 :/* 13.2 */ 0962 add_triangle( tiling13_2_[i_config][0], 6 ) ; 0963 break ; 0964 case 40 :/* 13.2 */ 0965 add_triangle( tiling13_2_[i_config][1], 6 ) ; 0966 break ; 0967 case 41 :/* 13.2 */ 0968 add_triangle( tiling13_2_[i_config][2], 6 ) ; 0969 break ; 0970 case 42 :/* 13.2 */ 0971 add_triangle( tiling13_2_[i_config][3], 6 ) ; 0972 break ; 0973 case 43 :/* 13.2 */ 0974 add_triangle( tiling13_2_[i_config][4], 6 ) ; 0975 break ; 0976 case 44 :/* 13.2 */ 0977 add_triangle( tiling13_2_[i_config][5], 6 ) ; 0978 break ; 0979 0980 case 45 :/* 13.1 */ 0981 add_triangle( tiling13_1_[i_config], 4 ) ; 0982 break ; 0983 0984 default : 0985 // Aqsis::log() << warning << "Impossible case 13 ?" << std::endl; 0986 print_cube() ; 0987 } 0988 break ; 0989 0990 case 14 : 0991 add_triangle( tiling14[i_config], 4 ) ; 0992 break ; 0993 }; 0994 } 0995 //_____________________________________________________________________________ 0996 0997 0998 0999 //_____________________________________________________________________________ 1000 // Adding triangles 1001 void MarchingCubes::add_triangle( const int* trig, char n, int v12 ) 1002 //----------------------------------------------------------------------------- 1003 { 1004 int tv[3] ; 1005 1006 for( int t = 0 ; t < 3*n ; t++ ) 1007 { 1008 int t3 = t % 3; 1009 switch( trig[t] ) 1010 { 1011 case 0 : 1012 tv[ t3 ] = get_x_vert( i_i , i_j , i_k ) ; 1013 break ; 1014 case 1 : 1015 tv[ t3 ] = get_y_vert(i_i+1, i_j , i_k ) ; 1016 break ; 1017 case 2 : 1018 tv[ t3 ] = get_x_vert( i_i ,i_j+1, i_k ) ; 1019 break ; 1020 case 3 : 1021 tv[ t3 ] = get_y_vert( i_i , i_j , i_k ) ; 1022 break ; 1023 case 4 : 1024 tv[ t3 ] = get_x_vert( i_i , i_j ,i_k+1) ; 1025 break ; 1026 case 5 : 1027 tv[ t3 ] = get_y_vert(i_i+1, i_j ,i_k+1) ; 1028 break ; 1029 case 6 : 1030 tv[ t3 ] = get_x_vert( i_i ,i_j+1,i_k+1) ; 1031 break ; 1032 case 7 : 1033 tv[ t3 ] = get_y_vert( i_i , i_j ,i_k+1) ; 1034 break ; 1035 case 8 : 1036 tv[ t3 ] = get_z_vert( i_i , i_j , i_k ) ; 1037 break ; 1038 case 9 : 1039 tv[ t3 ] = get_z_vert(i_i+1, i_j , i_k ) ; 1040 break ; 1041 case 10 : 1042 tv[ t3 ] = get_z_vert(i_i+1,i_j+1, i_k ) ; 1043 break ; 1044 case 11 : 1045 tv[ t3 ] = get_z_vert( i_i ,i_j+1, i_k ) ; 1046 break ; 1047 case 12 : 1048 tv[ t3 ] = v12 ; 1049 break ; 1050 default : 1051 break ; 1052 } 1053 1054 if( tv[t3] == -1 ) 1055 { 1056 // Aqsis::log() << warning << "Invalid triangle " << i_ntrigs << std::endl; 1057 print_cube() ; 1058 } 1059 1060 if( t3 == 2 ) 1061 { 1062 if( i_ntrigs >= i_Ntrigs ) 1063 { 1064 Triangle *temp = i_triangles ; 1065 i_triangles = new Triangle[ i_ntrigs + 1024] ; 1066 memcpy( i_triangles, temp, i_Ntrigs*sizeof(Triangle) ) ; 1067 delete[] temp ; 1068 /* 1069 Aqsis::log() << warning << "allocated triangles " << i_Ntrigs << std::endl; 1070 */ 1071 i_Ntrigs = i_ntrigs + 1024 ; 1072 } 1073 1074 Triangle *T = i_triangles + i_ntrigs++ ; 1075 T->v1 = tv[0] ; 1076 T->v2 = tv[1] ; 1077 T->v3 = tv[2] ; 1078 } 1079 } 1080 } 1081 //_____________________________________________________________________________ 1082 1083 1084 1085 //_____________________________________________________________________________ 1086 // Calculating gradient 1087 1088 double MarchingCubes::get_x_grad( const int i, const int j, const int k ) 1089 //----------------------------------------------------------------------------- 1090 { 1091 if( i > 0 ) 1092 { 1093 if ( i < i_size_x - 1 ) 1094 return ( get_data( i+1, j, k ) - get_data( i-1, j, k ) ) / 2 ; 1095 else 1096 return get_data( i, j, k ) - get_data( i-1, j, k ) ; 1097 } 1098 else 1099 return get_data( i+1, j, k ) - get_data( i, j, k ) ; 1100 } 1101 //----------------------------------------------------------------------------- 1102 1103 double MarchingCubes::get_y_grad( const int i, const int j, const int k ) 1104 //----------------------------------------------------------------------------- 1105 { 1106 if( j > 0 ) 1107 { 1108 if ( j < i_size_y - 1 ) 1109 return ( get_data( i, j+1, k ) - get_data( i, j-1, k ) ) / 2 ; 1110 else 1111 return get_data( i, j, k ) - get_data( i, j-1, k ) ; 1112 } 1113 else 1114 return get_data( i, j+1, k ) - get_data( i, j, k ) ; 1115 } 1116 //----------------------------------------------------------------------------- 1117 1118 double MarchingCubes::get_z_grad( const int i, const int j, const int k ) 1119 //----------------------------------------------------------------------------- 1120 { 1121 if( k > 0 ) 1122 { 1123 if ( k < i_size_z - 1 ) 1124 return ( get_data( i, j, k+1 ) - get_data( i, j, k-1 ) ) / 2 ; 1125 else 1126 return get_data( i, j, k ) - get_data( i, j, k-1 ) ; 1127 } 1128 else 1129 return get_data( i, j, k+1 ) - get_data( i, j, k ) ; 1130 } 1131 //_____________________________________________________________________________ 1132 1133 1134 //_____________________________________________________________________________ 1135 // Adding vertices 1136 1137 void MarchingCubes::test_vertex_addition() 1138 { 1139 if( i_nverts >= i_Nverts ) 1140 { 1141 Vertex *temp = i_vertices ; 1142 i_vertices = new Vertex[ i_nverts + 1024] ; 1143 memcpy( i_vertices, temp, i_Nverts*sizeof(Vertex) ) ; 1144 delete[] temp ; 1145 /* 1146 Aqsis::log() << warning << "allocated vertices " << i_Nverts << std::endl; 1147 */ 1148 i_Nverts = i_nverts + 1024 ; 1149 } 1150 } 1151 1152 1153 int MarchingCubes::add_x_vertex( ) 1154 //----------------------------------------------------------------------------- 1155 { 1156 test_vertex_addition() ; 1157 Vertex *vert = i_vertices + i_nverts++ ; 1158 1159 double u = ( i_cube[0] ) / ( i_cube[0] - i_cube[1] ) ; 1160 1161 //el mundo de analitza no es discreto, por eso los comentarios al cod original 1162 //una cosa mas --- el u es el resultado de la interpolacion lineal sobre el edge 1163 // vert->x = (float)i_i+u; 1164 // vert->y = (float) i_j ; 1165 // vert->z = (float) i_k ; 1166 1167 vert->x = xmin+hx*(i_i+u); 1168 vert->y = ymin+hy*i_j; 1169 vert->z = zmin+hz*i_k; 1170 1171 1172 #ifdef COMPUTE_NORMALS 1173 vert->nx = (1-u)*get_x_grad(i_i,i_j,i_k) + u*get_x_grad(i_i+1,i_j,i_k) ; 1174 vert->ny = (1-u)*get_y_grad(i_i,i_j,i_k) + u*get_y_grad(i_i+1,i_j,i_k) ; 1175 vert->nz = (1-u)*get_z_grad(i_i,i_j,i_k) + u*get_z_grad(i_i+1,i_j,i_k) ; 1176 1177 u = (float) sqrt( vert->nx * vert->nx + vert->ny * vert->ny +vert->nz * vert->nz ) ; 1178 if( u > 0 ) 1179 { 1180 vert->nx /= u ; 1181 vert->ny /= u ; 1182 vert->nz /= u ; 1183 } 1184 #endif 1185 1186 return i_nverts-1 ; 1187 } 1188 //----------------------------------------------------------------------------- 1189 1190 int MarchingCubes::add_y_vertex( ) 1191 //----------------------------------------------------------------------------- 1192 { 1193 test_vertex_addition() ; 1194 Vertex *vert = i_vertices + i_nverts++ ; 1195 1196 double u = ( i_cube[0] ) / ( i_cube[0] - i_cube[3] ) ; 1197 1198 // vert->x = (float) i_i ; 1199 // vert->y = (float)i_j+u; 1200 // vert->z = (float) i_k ; 1201 1202 vert->x = xmin+hx*i_i; 1203 vert->y = ymin+hy*(i_j+u); 1204 vert->z = zmin+hz*i_k; 1205 1206 1207 1208 #ifdef COMPUTE_NORMALS 1209 vert->nx = (1-u)*get_x_grad(i_i,i_j,i_k) + u*get_x_grad(i_i,i_j+1,i_k) ; 1210 vert->ny = (1-u)*get_y_grad(i_i,i_j,i_k) + u*get_y_grad(i_i,i_j+1,i_k) ; 1211 vert->nz = (1-u)*get_z_grad(i_i,i_j,i_k) + u*get_z_grad(i_i,i_j+1,i_k) ; 1212 1213 u = (float) sqrt( vert->nx * vert->nx + vert->ny * vert->ny +vert->nz * vert->nz ) ; 1214 if( u > 0 ) 1215 { 1216 vert->nx /= u ; 1217 vert->ny /= u ; 1218 vert->nz /= u ; 1219 } 1220 #endif 1221 1222 return i_nverts-1 ; 1223 } 1224 //----------------------------------------------------------------------------- 1225 1226 int MarchingCubes::add_z_vertex( ) 1227 //----------------------------------------------------------------------------- 1228 { 1229 test_vertex_addition() ; 1230 Vertex *vert = i_vertices + i_nverts++ ; 1231 1232 double u = ( i_cube[0] ) / ( i_cube[0] - i_cube[4] ) ; 1233 1234 // vert->x = (float) i_i ; 1235 // vert->y = (float) i_j ; 1236 // vert->z = (float)i_k+u; 1237 1238 vert->x = xmin+hx*i_i; 1239 vert->y = ymin+hy*i_j; 1240 vert->z = zmin+hz*(i_k+u); 1241 1242 1243 #ifdef COMPUTE_NORMALS 1244 vert->nx = (1-u)*get_x_grad(i_i,i_j,i_k) + u*get_x_grad(i_i,i_j,i_k+1) ; 1245 vert->ny = (1-u)*get_y_grad(i_i,i_j,i_k) + u*get_y_grad(i_i,i_j,i_k+1) ; 1246 vert->nz = (1-u)*get_z_grad(i_i,i_j,i_k) + u*get_z_grad(i_i,i_j,i_k+1) ; 1247 1248 u = (float) sqrt( vert->nx * vert->nx + vert->ny * vert->ny +vert->nz * vert->nz ) ; 1249 if( u > 0 ) 1250 { 1251 vert->nx /= u ; 1252 vert->ny /= u ; 1253 vert->nz /= u ; 1254 } 1255 #endif 1256 1257 return i_nverts-1 ; 1258 } 1259 1260 1261 int MarchingCubes::add_c_vertex( ) 1262 //----------------------------------------------------------------------------- 1263 { 1264 test_vertex_addition() ; 1265 Vertex *vert = i_vertices + i_nverts++ ; 1266 1267 double u = 0 ; 1268 int vid ; 1269 1270 vert->x = vert->y = vert->z = 0; 1271 #ifdef COMPUTE_NORMALS 1272 vert->nx = vert->ny = vert->nz = 0 ; 1273 #endif 1274 1275 // Computes the average of the intersection points of the cube 1276 vid = get_x_vert( i_i , i_j , i_k ) ; 1277 if( vid != -1 ) 1278 { 1279 ++u ; 1280 const Vertex &v = i_vertices[vid] ; 1281 vert->x += v.x ; 1282 vert->y += v.y ; 1283 vert->z += v.z ; 1284 #ifdef COMPUTE_NORMALS 1285 vert->nx += v.nx ; 1286 vert->ny += v.ny ; 1287 vert->nz += v.nz ; 1288 #endif 1289 } 1290 vid = get_y_vert(i_i+1, i_j , i_k ) ; 1291 if( vid != -1 ) 1292 { 1293 ++u ; 1294 const Vertex &v = i_vertices[vid] ; 1295 vert->x += v.x ; 1296 vert->y += v.y ; 1297 vert->z += v.z ; 1298 #ifdef COMPUTE_NORMALS 1299 vert->nx += v.nx ; 1300 vert->ny += v.ny ; 1301 vert->nz += v.nz ; 1302 #endif 1303 } 1304 vid = get_x_vert( i_i ,i_j+1, i_k ) ; 1305 if( vid != -1 ) 1306 { 1307 ++u ; 1308 const Vertex &v = i_vertices[vid] ; 1309 vert->x += v.x ; 1310 vert->y += v.y ; 1311 vert->z += v.z ; 1312 #ifdef COMPUTE_NORMALS 1313 vert->nx += v.nx ; 1314 vert->ny += v.ny ; 1315 vert->nz += v.nz ; 1316 #endif 1317 } 1318 vid = get_y_vert( i_i , i_j , i_k ) ; 1319 if( vid != -1 ) 1320 { 1321 ++u ; 1322 const Vertex &v = i_vertices[vid] ; 1323 vert->x += v.x ; 1324 vert->y += v.y ; 1325 vert->z += v.z ; 1326 #ifdef COMPUTE_NORMALS 1327 vert->nx += v.nx ; 1328 vert->ny += v.ny ; 1329 vert->nz += v.nz ; 1330 #endif 1331 } 1332 vid = get_x_vert( i_i , i_j ,i_k+1) ; 1333 if( vid != -1 ) 1334 { 1335 ++u ; 1336 const Vertex &v = i_vertices[vid] ; 1337 vert->x += v.x ; 1338 vert->y += v.y ; 1339 vert->z += v.z ; 1340 #ifdef COMPUTE_NORMALS 1341 vert->nx += v.nx ; 1342 vert->ny += v.ny ; 1343 vert->nz += v.nz ; 1344 #endif 1345 } 1346 vid = get_y_vert(i_i+1, i_j ,i_k+1) ; 1347 if( vid != -1 ) 1348 { 1349 ++u ; 1350 const Vertex &v = i_vertices[vid] ; 1351 vert->x += v.x ; 1352 vert->y += v.y ; 1353 vert->z += v.z ; 1354 #ifdef COMPUTE_NORMALS 1355 vert->nx += v.nx ; 1356 vert->ny += v.ny ; 1357 vert->nz += v.nz ; 1358 #endif 1359 } 1360 vid = get_x_vert( i_i ,i_j+1,i_k+1) ; 1361 if( vid != -1 ) 1362 { 1363 ++u ; 1364 const Vertex &v = i_vertices[vid] ; 1365 vert->x += v.x ; 1366 vert->y += v.y ; 1367 vert->z += v.z ; 1368 #ifdef COMPUTE_NORMALS 1369 vert->nx += v.nx ; 1370 vert->ny += v.ny ; 1371 vert->nz += v.nz ; 1372 #endif 1373 } 1374 vid = get_y_vert( i_i , i_j ,i_k+1) ; 1375 if( vid != -1 ) 1376 { 1377 ++u ; 1378 const Vertex &v = i_vertices[vid] ; 1379 vert->x += v.x ; 1380 vert->y += v.y ; 1381 vert->z += v.z ; 1382 #ifdef COMPUTE_NORMALS 1383 vert->nx += v.nx ; 1384 vert->ny += v.ny ; 1385 vert->nz += v.nz ; 1386 #endif 1387 } 1388 vid = get_z_vert( i_i , i_j , i_k ) ; 1389 if( vid != -1 ) 1390 { 1391 ++u ; 1392 const Vertex &v = i_vertices[vid] ; 1393 vert->x += v.x ; 1394 vert->y += v.y ; 1395 vert->z += v.z ; 1396 #ifdef COMPUTE_NORMALS 1397 vert->nx += v.nx ; 1398 vert->ny += v.ny ; 1399 vert->nz += v.nz ; 1400 #endif 1401 } 1402 vid = get_z_vert(i_i+1, i_j , i_k ) ; 1403 if( vid != -1 ) 1404 { 1405 ++u ; 1406 const Vertex &v = i_vertices[vid] ; 1407 vert->x += v.x ; 1408 vert->y += v.y ; 1409 vert->z += v.z ; 1410 #ifdef COMPUTE_NORMALS 1411 vert->nx += v.nx ; 1412 vert->ny += v.ny ; 1413 vert->nz += v.nz ; 1414 #endif 1415 } 1416 vid = get_z_vert(i_i+1,i_j+1, i_k ) ; 1417 if( vid != -1 ) 1418 { 1419 ++u ; 1420 const Vertex &v = i_vertices[vid] ; 1421 vert->x += v.x ; 1422 vert->y += v.y ; 1423 vert->z += v.z ; 1424 #ifdef COMPUTE_NORMALS 1425 vert->nx += v.nx ; 1426 vert->ny += v.ny ; 1427 vert->nz += v.nz ; 1428 #endif 1429 } 1430 vid = get_z_vert( i_i ,i_j+1, i_k ) ; 1431 if( vid != -1 ) 1432 { 1433 ++u ; 1434 const Vertex &v = i_vertices[vid] ; 1435 vert->x += v.x ; 1436 vert->y += v.y ; 1437 vert->z += v.z ; 1438 #ifdef COMPUTE_NORMALS 1439 vert->nx += v.nx ; 1440 vert->ny += v.ny ; 1441 vert->nz += v.nz ; 1442 #endif 1443 } 1444 1445 vert->x /= u ; 1446 vert->y /= u ; 1447 vert->z /= u ; 1448 1449 #ifdef COMPUTE_NORMALS 1450 u = (float) sqrt( vert->nx * vert->nx + vert->ny * vert->ny +vert->nz * vert->nz ) ; 1451 if( u > 0 ) 1452 { 1453 vert->nx /= u ; 1454 vert->ny /= u ; 1455 vert->nz /= u ; 1456 } 1457 #endif 1458 1459 return i_nverts-1 ; 1460 } 1461 //_____________________________________________________________________________ 1462 1463 1464 1465 //_____________________________________________________________________________ 1466 //_____________________________________________________________________________ 1467 1468 1469 1470 1471 void MarchingCubes::write(const char *fn, bool /*bin */) 1472 //----------------------------------------------------------------------------- 1473 { 1474 FILE *fp = fopen( fn, "w" ); 1475 fprintf(fp, "%d %d\n", i_nverts, i_ntrigs); 1476 1477 int i ; 1478 1479 for ( i = 0; i < i_nverts; i++ ) 1480 fprintf(fp, "%f %f %f\n", i_vertices[i].x, i_vertices[i].y, i_vertices[i].z); 1481 1482 for ( i = 0; i < i_ntrigs; i++ ) 1483 { 1484 fprintf(fp, "%d %d %d \n", i_triangles[i].v1, i_triangles[i].v2, i_triangles[i].v3); 1485 } 1486 1487 fclose( fp ) ; 1488 } 1489 //_____________________________________________________________________________ 1490