00001 #include "Propagation.h" 00002 00003 //------------------------------------------------------------------------------------------------------------------------------------------ 00004 //CLASS: Vector ------------------------------------------------------------------------------------------------------------------------ 00005 //------------------------------------------------------------------------------------------------------------------------------------------ 00006 //Constructor 00007 Vector::Vector() 00008 { 00009 _plane = -1; 00010 _var = -1; 00011 } 00012 //Destructor 00013 Vector::~Vector() 00014 { 00015 _vec.clear(); 00016 _vecX.clear(); 00017 _vecY.clear(); 00018 _vecZ.clear(); 00019 } 00020 00021 //------------------------------------------------------------------------------------------------------------------------------------------ 00022 void Vector::set_vec(double val) 00023 { 00024 _vec.push_back(val); 00025 } 00026 //------------------------------------------------------------------------------------------------------------------------------------------ 00027 void Vector::set_var(double val) 00028 { 00029 _var = val; 00030 } 00031 //------------------------------------------------------------------------------------------------------------------------------------------ 00032 double Vector::get_vec(int id) 00033 { 00034 if(_vec.size() != 0) 00035 { 00036 return _vec[id]; 00037 } 00038 return -1; 00039 } 00040 //------------------------------------------------------------------------------------------------------------------------------------------ 00041 double Vector::get_var() 00042 { 00043 return _var; 00044 } 00045 //------------------------------------------------------------------------------------------------------------------------------------------ 00046 int Vector::getsize_vec() 00047 { 00048 return _vec.size(); 00049 } 00050 //------------------------------------------------------------------------------------------------------------------------------------------ 00051 void Vector::copyVector( std::vector<Vector>*vec1, std::vector<Vector>*vec2 ) 00052 { 00053 int i,j; 00054 if( vec1->size() != 0 ) 00055 { 00056 vec2->clear(); 00057 for(i=0; i<(int)(vec1->size()); i++) 00058 { 00059 Vector *temp = new Vector(); 00060 temp->set_var( (*vec1)[i].get_var() ); 00061 for(j=0; j<(*vec1)[i].getsize_vec(); j++) 00062 { 00063 temp->set_vec( (*vec1)[i].get_vec(j) ); 00064 } 00065 vec2->push_back(*temp); 00066 delete temp; 00067 } 00068 } 00069 } 00070 //------------------------------------------------------------------------------------------------------------------------------------------ 00071 void Vector::printVector(std::vector<Vector>*vec1) 00072 { 00073 int i,j; 00074 for(i=0; i<(int)(vec1->size()); i++) 00075 { 00076 printf("\n Pos (%d) => var = %f",i,(*vec1)[i].get_var()); 00077 for(j=0; j<(*vec1)[i].getsize_vec(); j++) 00078 { 00079 printf("\n vec(%d) = %f",j,(*vec1)[i].get_vec(j)); 00080 } 00081 } 00082 } 00083 //------------------------------------------------------------------------------------------------------------------------------------------ 00084 void Vector::set_plane(int val) 00085 { 00086 _plane = val; 00087 } 00088 //------------------------------------------------------------------------------------------------------------------------------------------ 00089 void Vector::set_x(double val) 00090 { 00091 _vecX.push_back(val); 00092 } 00093 //------------------------------------------------------------------------------------------------------------------------------------------ 00094 //------------------------------------------------------------------------------------------------------------------------------------------ 00095 void Vector::set_y(double val) 00096 { 00097 _vecY.push_back(val); 00098 } 00099 //------------------------------------------------------------------------------------------------------------------------------------------ 00100 //------------------------------------------------------------------------------------------------------------------------------------------ 00101 void Vector::set_z(double val) 00102 { 00103 _vecZ.push_back(val); 00104 } 00105 //------------------------------------------------------------------------------------------------------------------------------------------ 00106 int Vector::get_plane() 00107 { 00108 return _plane; 00109 } 00110 //------------------------------------------------------------------------------------------------------------------------------------------ 00111 double Vector::get_x(int id) 00112 { 00113 if( (-1<id) && (id<(int)(_vecX.size()) ) ) 00114 { 00115 return _vecX[id]; 00116 } 00117 return -1; 00118 } 00119 //------------------------------------------------------------------------------------------------------------------------------------------ 00120 int Vector::getsize_x() 00121 { 00122 if(_vecX.size() != 0) 00123 { 00124 return _vecX.size(); 00125 } 00126 return -1; 00127 } 00128 //------------------------------------------------------------------------------------------------------------------------------------------ 00129 double Vector::get_y(int id) 00130 { 00131 if( (-1<id) && (id<(int)(_vecY.size())) ) 00132 { 00133 return _vecY[id]; 00134 } 00135 return -1; 00136 } 00137 //------------------------------------------------------------------------------------------------------------------------------------------ 00138 int Vector::getsize_y() 00139 { 00140 if(_vecY.size() != 0) 00141 { 00142 return _vecY.size(); 00143 } 00144 return -1; 00145 } 00146 //------------------------------------------------------------------------------------------------------------------------------------------ 00147 int Vector::getsize_z() 00148 { 00149 if(_vecZ.size() != 0) 00150 { 00151 return _vecZ.size(); 00152 } 00153 return -1; 00154 } 00155 //------------------------------------------------------------------------------------------------------------------------------------------ 00156 double Vector::get_z(int id) 00157 { 00158 if( (-1<id) && (id<(int)(_vecZ.size())) ) 00159 { 00160 return _vecZ[id]; 00161 } 00162 return -1; 00163 } 00164 //------------------------------------------------------------------------------------------------------------------------------------------ 00165 std::vector<double> Vector::getVec() 00166 { 00167 return _vec; 00168 } 00169 //------------------------------------------------------------------------------------------------------------------------------------------ 00170 void Vector::resetVec() 00171 { 00172 _vec.clear(); 00173 } 00174 //------------------------------------------------------------------------------------------------------------------------------------------ 00175 //------------------------------------------------------------------------------------------------------------------------------------------ 00176 //------------------------------------------------------------------------------------------------------------------------------------------ 00177 00178 00179 //Constructor 00180 PropContour::PropContour() 00181 { 00182 _interpnumber = 100; 00183 } 00184 00185 //Destructor 00186 PropContour::~PropContour() 00187 { 00188 ResetPlaneVector(); 00189 ResetKeyContours(); 00190 } 00191 00192 /* EED 03/07/2008 00193 //------------------------------------------------------------------------------------ 00194 double PropContour::RBF_WendLand(double norm, double m_rad) 00195 { 00196 double y; 00197 norm = norm / m_rad; 00198 y = pow( 1-norm,4 ) * ( (4*norm) + 1 ); 00199 if(norm >= 1) 00200 { 00201 y = 0; 00202 } 00203 return y; 00204 } 00205 //------------------------------------------------------------------------------------ 00206 double PropContour::RBF_ThinPlate(double norm) 00207 { 00208 double y; 00209 if(norm == 0) 00210 { 00211 y = 0; 00212 } 00213 else 00214 { 00215 y = pow(norm,2)*log(norm); 00216 } 00217 return y; 00218 } 00219 //------------------------------------------------------------------------------------ 00220 vtkImageData* PropContour::method_RBF ( double rad, std::vector<double>*CoordX, std::vector<double>*CoordY, 00221 std::vector<double>*CoordZ ) 00222 { 00223 _dimImage[0] = 100; // X Axis 00224 _dimImage[1] = 100; // Y Axis 00225 _dimImage[2] = 1; // Z Axis 00226 int pointsSize = CoordX->size(); 00227 double spc[3]={1,1,1}; 00228 double norm = 0, val = 0; 00229 00230 vnl_matrix<double> A( pointsSize,pointsSize, 0.0 ); 00231 vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 ); 00232 vnl_matrix<double> I( pointsSize,pointsSize, 0.0 ); 00233 vnl_vector<double> H( pointsSize, 1 ); 00234 vnl_vector<double> D( pointsSize, 0.0 ); 00235 vnl_vector<double> Q( 2, 0.0 ); 00236 vnl_matrix<double> Impl( _dimImage[0],_dimImage[1], 0.0 ); 00237 00238 unsigned short *pValue; 00239 imagedataValue = vtkImageData::New(); 00240 imagedataValue->SetScalarTypeToUnsignedShort(); 00241 imagedataValue->SetSpacing(spc); 00242 imagedataValue->SetDimensions(_dimImage); 00243 imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1); 00244 imagedataValue->AllocateScalars(); 00245 imagedataValue->Update(); 00246 00247 int i,j,h; 00248 for(i=0; i<pointsSize; i++) 00249 { 00250 for(j=0; j<pointsSize; j++) 00251 { 00252 norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) ); 00253 A(i,j) = RBF_WendLand(norm,rad); 00254 //A(i,j) = RBF_ThinPlate(norm); 00255 } 00256 } 00257 00258 Ainv = vnl_matrix_inverse<double>(A); 00259 D = Ainv* H; 00260 00261 for(i=0; i < _dimImage[1]; i++) 00262 { 00263 for(j=0; j < _dimImage[0]; j++) 00264 { 00265 Q(0) = i; 00266 Q(1) = j; 00267 val = 0; 00268 for(h=0; h<pointsSize; h++) 00269 { 00270 norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) ); 00271 //val = val + ( D(h) * RBF_ThinPlate(norm) ); 00272 val = val + ( D(h) * RBF_WendLand(norm,rad) ); 00273 } 00274 Impl(i,j) = val - 1.0; 00275 00276 // if ( (Impl(i,j)>=-0.001) && (Impl(i,j) <=0.001) ) 00277 // {Impl(i,j)=128; 00278 // } 00279 // else 00280 // {Impl(i,j)=0;} 00281 00282 00283 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,0); 00284 // vvalue = (Impl(i,j)*256+256); 00285 // if (vvalue) < 0) 00286 // { 00287 // *pValue = 0; 00288 // } else { 00289 // *pValue = vvalue; 00290 // } 00291 00292 pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,0); 00293 if ( Impl(i,j) >=0 ) 00294 { 00295 *pValue = 128; 00296 } 00297 else 00298 { 00299 *pValue = 0; 00300 } 00301 } 00302 } 00303 00304 for(i=0; i<pointsSize; i=i+5) 00305 { 00306 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,0); 00307 *pValue=255; 00308 } 00309 00310 return imagedataValue; 00311 } 00312 //---------------------------------------------------------------------------------------------------- 00313 double PropContour::RBF_ThinPlate_3D(double norm) 00314 { 00315 return norm = pow( norm,3 ); 00316 } 00317 //---------------------------------------------------------------------------------------------------- 00318 vtkImageData* PropContour::method_RBF_3D ( double rad, std::vector<double>*CoordX, std::vector<double>*CoordY, 00319 std::vector<double>*CoordZ ) 00320 { 00321 long interval = wxGetElapsedTime(TRUE); 00322 00323 int i,j,k,h; 00324 int pointsSize = CoordX->size(); 00325 double max = -1, min = 1000, minz = 100000, maxz = -10000; 00326 00327 for( i=0; i<pointsSize; i++ ) 00328 { 00329 if( (*CoordX)[i] > max ) 00330 { 00331 max = (*CoordX)[i]; 00332 } 00333 if( (*CoordY)[i] > max ) 00334 { 00335 max = (*CoordY)[i]; 00336 } 00337 if( (*CoordX)[i] < min ) 00338 { 00339 min = (*CoordX)[i]; 00340 } 00341 if( (*CoordY)[i] < min ) 00342 { 00343 min = (*CoordY)[i]; 00344 } 00345 if( (*CoordZ)[i] < minz ) 00346 { 00347 minz = (*CoordZ)[i]; 00348 } 00349 if( (*CoordZ)[i] > maxz ) 00350 { 00351 maxz = (*CoordZ)[i]; 00352 } 00353 } 00354 00355 _dimImage[0] = 200; // X axis 00356 _dimImage[1] = 200; // Y axis 00357 _dimImage[2] = 200; // Z axis 00358 00359 double spc[3]={1,1,1}; 00360 double norm = 0, val, vvalue; 00361 00362 vnl_matrix<double> A( pointsSize,pointsSize, 0.0 ); 00363 vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 ); 00364 vnl_matrix<double> I( pointsSize,pointsSize, 0.0 ); 00365 vnl_vector<double> H( pointsSize, 1 ); 00366 vnl_vector<double> D( pointsSize, 0.0 ); 00367 vnl_vector<double> Q( 3, 0.0 ); 00368 00369 unsigned short *pValue; 00370 imagedataValue = vtkImageData::New(); 00371 imagedataValue->SetScalarTypeToUnsignedShort(); 00372 imagedataValue->SetSpacing(spc); 00373 imagedataValue->SetDimensions(_dimImage); 00374 imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1); 00375 imagedataValue->AllocateScalars(); 00376 imagedataValue->Update(); 00377 00378 for(i=0; i<pointsSize; i++) 00379 { 00380 for(j=0; j<pointsSize; j++) 00381 { 00382 norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) + pow((*CoordZ)[i]-(*CoordZ)[j],2) ); 00383 A(i,j) = RBF_WendLand(norm,rad); 00384 } 00385 } 00386 D = vnl_matrix_inverse<double>(A)* H; 00387 00388 //Inicialization 00389 // for(k=0; k<_dimImage[2]; k++) 00390 // { 00391 // for(j=0; j<_dimImage[1]; j++) 00392 // { 00393 // for(i=0; i<_dimImage[0]; i++) 00394 // { 00395 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k); 00396 // *pValue = 0; 00397 // } 00398 // } 00399 // } 00400 00401 //Filling 00402 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(0,0,0); 00403 // i = 0; 00404 // j = 0; 00405 // k = 0; 00406 // for(h=0; h<pointsSize; h++) 00407 // { 00408 // norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 )); 00409 // //val = val + ( D(h) * RBF_ThinPlate_3D(norm) ); 00410 // val = val + ( D(h) * RBF_WendLand(norm,rad) ); 00411 // if( h == pointsSize-1 ) 00412 // { 00413 // val = val - 1; 00414 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k); 00415 // if ( val >=0 ) 00416 // { 00417 // *pValue = 128; 00418 // } 00419 // else 00420 // { 00421 // *pValue = 0; 00422 // } 00423 // //pValue++; 00424 // j++; 00425 // h = -1; 00426 // val = 0; 00427 // if( j == _dimImage[0] ) 00428 // { 00429 // i++; 00430 // j = 0; 00431 // h = -1; 00432 // if( i == _dimImage[1]) 00433 // { 00434 // k++; 00435 // i = 0; 00436 // j = 0; 00437 // h = -1; 00438 // if(k == _dimImage[2]) 00439 // { 00440 // h = pointsSize; 00441 // } 00442 // } 00443 // } 00444 // } 00445 // } 00446 00447 i = (int)min-10; 00448 j = (int)min-10; 00449 k = (int)minz-10; 00450 val = 0; 00451 00452 for(h=0; h<pointsSize; h++) 00453 { 00454 norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 )); 00455 val = val + ( D(h) * RBF_WendLand(norm,rad) ); 00456 if( h == pointsSize-1 ) 00457 { 00458 val = val - 1; 00459 pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k); 00460 if ( val >=0 ) 00461 { 00462 *pValue = 128; 00463 } 00464 else 00465 { 00466 *pValue = 0; 00467 } 00468 //pValue++; 00469 j++; 00470 h = -1; 00471 val = 0; 00472 if( j == (int)max+10) 00473 { 00474 i++; 00475 j = min-10; 00476 h = -1; 00477 if( i == (int)max+10) 00478 { 00479 k++; 00480 i = min-10; 00481 j = min-10; 00482 h = -1; 00483 if(k == maxz+10) 00484 { 00485 h = pointsSize; 00486 } 00487 } 00488 } 00489 } 00490 } 00491 00492 interval = wxGetElapsedTime(FALSE); 00493 long interPlane = interval/_dimImage[2]; 00494 00495 for(i=0; i<pointsSize; i++) 00496 { 00497 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,(int)(*CoordZ)[i]); 00498 *pValue=255; 00499 } 00500 00501 long intervalPC = wxGetElapsedTime(); 00502 00503 printf("\n\n JSTG - PropContour::method_RBF_3D ------------------------"); 00504 printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %lld (ms)",interval); 00505 printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC); 00506 printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %lld (ms)",interPlane); 00507 printf("\n NUMBER OF PLANES................................. %d",k); 00508 printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize); 00509 printf("\n ------------------------------------------------------------"); 00510 00511 return imagedataValue; 00512 } 00513 //---------------------------------------------------------------------------------------------------- 00514 00515 vtkImageData* PropContour::method_RBF_3D_ThinPlate ( double rad, std::vector<double>*CoordX, std::vector<double>*CoordY, 00516 std::vector<double>*CoordZ ) 00517 { 00518 long interval = wxGetElapsedTime(TRUE); 00519 00520 int i,j,k,h; 00521 int pointsSize = CoordX->size(); 00522 double max = -1, min = 1000, minz = 100000, maxz = -10000; 00523 00524 for( i=0; i<pointsSize; i++ ) 00525 { 00526 if( (*CoordX)[i] > max ) 00527 { 00528 max = (*CoordX)[i]; 00529 } 00530 if( (*CoordY)[i] > max ) 00531 { 00532 max = (*CoordY)[i]; 00533 } 00534 if( (*CoordX)[i] < min ) 00535 { 00536 min = (*CoordX)[i]; 00537 } 00538 if( (*CoordY)[i] < min ) 00539 { 00540 min = (*CoordY)[i]; 00541 } 00542 if( (*CoordZ)[i] < minz ) 00543 { 00544 minz = (*CoordZ)[i]; 00545 } 00546 if( (*CoordZ)[i] > maxz ) 00547 { 00548 maxz = (*CoordZ)[i]; 00549 } 00550 } 00551 00552 _dimImage[0] = 190; // X axis 00553 _dimImage[1] = 190; // Y axis 00554 _dimImage[2] = 190; // Z axis 00555 00556 double spc[3]={1,1,1}; 00557 double norm = 0, val, vvalue; 00558 00559 vnl_matrix<double> A( pointsSize,pointsSize, 0.0 ); 00560 vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 ); 00561 vnl_matrix<double> I( pointsSize,pointsSize, 0.0 ); 00562 vnl_vector<double> H( pointsSize, 1 ); 00563 vnl_vector<double> D( pointsSize, 0.0 ); 00564 vnl_vector<double> Q( 3, 0.0 ); 00565 00566 unsigned short *pValue; 00567 imagedataValue = vtkImageData::New(); 00568 imagedataValue->SetScalarTypeToUnsignedShort(); 00569 imagedataValue->SetSpacing(spc); 00570 imagedataValue->SetDimensions(_dimImage); 00571 imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,-10,_dimImage[2]-1); 00572 imagedataValue->AllocateScalars(); 00573 imagedataValue->Update(); 00574 00575 for(i=0; i<pointsSize; i++) 00576 { 00577 for(j=0; j<pointsSize; j++) 00578 { 00579 norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) + pow((*CoordZ)[i]-(*CoordZ)[j],2) ); 00580 A(i,j) = RBF_ThinPlate_3D(norm); 00581 } 00582 } 00583 D = vnl_matrix_inverse<double>(A)* H; 00584 00585 //Inicialization 00586 // for(k=0; k<_dimImage[2]; k++) 00587 // { 00588 // for(j=0; j<_dimImage[1]; j++) 00589 // { 00590 // for(i=0; i<_dimImage[0]; i++) 00591 // { 00592 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k); 00593 // *pValue = 0; 00594 // } 00595 // } 00596 // } 00597 00598 //Filling 00599 00600 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(0,0,0); 00601 // i = 0; 00602 // j = 0; 00603 // k = 0; 00604 // for(h=0; h<pointsSize; h++) 00605 // { 00606 // norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 )); 00607 // //val = val + ( D(h) * RBF_ThinPlate_3D(norm) ); 00608 // val = val + ( D(h) * RBF_WendLand(norm,rad) ); 00609 // if( h == pointsSize-1 ) 00610 // { 00611 // val = val - 1; 00612 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k); 00613 // if ( val >=0 ) 00614 // { 00615 // *pValue = 128; 00616 // } 00617 // else 00618 // { 00619 // *pValue = 0; 00620 // } 00621 // //pValue++; 00622 // j++; 00623 // h = -1; 00624 // val = 0; 00625 // if( j == _dimImage[0] ) 00626 // { 00627 // i++; 00628 // j = 0; 00629 // h = -1; 00630 // if( i == _dimImage[1]) 00631 // { 00632 // k++; 00633 // i = 0; 00634 // j = 0; 00635 // h = -1; 00636 // if(k == _dimImage[2]) 00637 // { 00638 // h = pointsSize; 00639 // } 00640 // } 00641 // } 00642 // } 00643 // } 00644 00645 i = (int)min-10; 00646 j = (int)min-10; 00647 k = (int)minz-10; 00648 val = 0; 00649 00650 for(h=0; h<pointsSize; h++) 00651 { 00652 norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 )); 00653 val = val + ( D(h) * RBF_ThinPlate_3D(norm) ); 00654 if( h == pointsSize-1 ) 00655 { 00656 val = val - 1; 00657 pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k); 00658 if ( val >=0 ) 00659 { 00660 *pValue = 128; 00661 } 00662 else 00663 { 00664 *pValue = 0; 00665 } 00666 //pValue++; 00667 j++; 00668 h = -1; 00669 val = 0; 00670 if( j == (int)max+10) 00671 { 00672 i++; 00673 j = min-10; 00674 h = -1; 00675 if( i == (int)max+10) 00676 { 00677 k++; 00678 i = min-10; 00679 j = min-10; 00680 h = -1; 00681 if(k == maxz+10) 00682 { 00683 h = pointsSize; 00684 } 00685 } 00686 } 00687 } 00688 } 00689 00690 interval = wxGetElapsedTime(FALSE); 00691 long interPlane = interval/_dimImage[2]; 00692 00693 for(i=0; i<pointsSize; i++) 00694 { 00695 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,(int)(*CoordZ)[i]); 00696 *pValue=255; 00697 } 00698 00699 long intervalPC = wxGetElapsedTime(); 00700 00701 printf("\n\n JSTG - PropContour::method_RBF_3D ------------------------"); 00702 printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %lld (ms)",interval); 00703 printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC); 00704 printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %lld (ms)",interPlane); 00705 printf("\n NUMBER OF PLANES................................. %d",k); 00706 printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize); 00707 printf("\n ------------------------------------------------------------"); 00708 00709 return imagedataValue; 00710 } 00711 00712 */ 00713 00714 //--------------------------------------------------------------------------------------------------------- 00715 void PropContour::ReadKeyContour(FILE* fd) 00716 { 00717 char firstline[30]; 00718 int size; 00719 double x,y; 00720 int z; 00721 std::vector<double> tempX; 00722 std::vector<double> tempY; 00723 std::vector<double> tempZ; 00724 tempX.clear(); 00725 tempY.clear(); 00726 tempZ.clear(); 00727 while(!feof(fd)) 00728 { 00729 //fscanf(fd," %s %d",&firstline,&size); // JPRx 00730 fscanf(fd," %s %d",firstline,&size); 00731 for(int i=0; i<size; i++) 00732 { 00733 fscanf(fd,"%lf %lf %d",&x,&y,&z); 00734 tempX.push_back(x); 00735 tempY.push_back(y); 00736 tempZ.push_back(z); 00737 } 00738 SetKeyContours(&tempX,&tempY,&tempZ); 00739 tempX.clear(); 00740 tempY.clear(); 00741 tempZ.clear(); 00742 } 00743 } 00744 //--------------------------------------------------------------------------------------------------------- 00745 int PropContour::VectorDirection(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ) 00746 { 00747 int dir,i; 00748 double SumX = 0,SumY = 0; 00749 double ax,ay,bx,by,axb; 00750 int size = InX->size(); 00751 for(i=0; i<size; i++) 00752 { 00753 SumX = SumX + (*InX)[i]; 00754 SumY = SumY + (*InY)[i]; 00755 } 00756 SumX = SumX/size; //Mass Center: X coord 00757 SumY = SumY/size; //Mass Center: Y coord 00758 00759 int positive = 0; 00760 int negative = 0; 00761 for(i=0; i<size; i++) 00762 { 00763 ax = (*InX)[i]-SumX; 00764 ay = (*InY)[i]-SumY; 00765 bx = (*InX)[(i+1)%size]-SumX; 00766 by = (*InY)[(i+1)%size]-SumY; 00767 axb = (ax*by) - (bx*ay); 00768 if(axb > 0) 00769 { 00770 positive++; 00771 } 00772 if(axb < 0) 00773 { 00774 negative++; 00775 } 00776 } 00777 if(positive >= negative) 00778 { 00779 dir = 1; 00780 } 00781 else 00782 { 00783 dir = -1; 00784 } 00785 00786 return dir; 00787 } 00788 //---------------------------------------------------------------------------------------------------- 00789 void PropContour::VectorOrder(int dir, int posini, std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ) 00790 { 00791 int size = InX->size(); 00792 int i; 00793 std::vector<double> tempX; 00794 std::vector<double> tempY; 00795 std::vector<double> tempZ; 00796 00797 tempX.clear(); 00798 tempY.clear(); 00799 tempZ.clear(); 00800 00801 for(i=0; i<size; i++) 00802 { 00803 if(dir == 1) 00804 { 00805 tempX.push_back((*InX)[posini]); 00806 tempY.push_back((*InY)[posini]); 00807 tempZ.push_back((*InZ)[posini]); 00808 posini++; 00809 if(posini == size) 00810 { 00811 posini = 0; 00812 } 00813 } 00814 if(dir == -1) 00815 { 00816 tempX.push_back((*InX)[posini]); 00817 tempY.push_back((*InY)[posini]); 00818 tempZ.push_back((*InZ)[posini]); 00819 posini--; 00820 if(posini < 0) 00821 { 00822 posini = size-1; 00823 } 00824 } 00825 } 00826 InX->clear(); 00827 InY->clear(); 00828 InZ->clear(); 00829 for(i=0; i<size; i++) 00830 { 00831 InX->push_back(tempX[i]); 00832 InY->push_back(tempY[i]); 00833 InZ->push_back(tempZ[i]); 00834 } 00835 00836 } 00837 //---------------------------------------------------------------------------------------------------- 00838 void PropContour::PreparePointsForSpline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ, 00839 std::vector<int>*Sizes) 00840 { 00841 int sizeS = Sizes->size(); 00842 //int sizeV = InX->size(); // JPRx 00843 int i,j,mem,posinic,dir,cont; 00844 double leX; 00845 00846 std::vector<double> tempX; 00847 std::vector<double> tempY; 00848 std::vector<double> tempZ; 00849 std::vector<double> lstX; 00850 std::vector<double> lstY; 00851 std::vector<double> lstZ; 00852 00853 lstX.clear(); 00854 lstY.clear(); 00855 lstZ.clear(); 00856 00857 mem = 0; 00858 cont = 0; 00859 for(i=0; i<sizeS; i++) 00860 { 00861 leX=1000; 00862 tempX.clear(); 00863 tempY.clear(); 00864 tempZ.clear(); 00865 for(j=0; j<(*Sizes)[i]; j++) 00866 { 00867 tempX.push_back((*InX)[j+mem]); 00868 tempY.push_back((*InY)[j+mem]); 00869 tempZ.push_back((*InZ)[j+mem]); 00870 if( (*InX)[j] < leX ) 00871 { 00872 posinic = j; 00873 leX = (*InX)[j]; 00874 } 00875 } 00876 mem = mem + (*Sizes)[i]; 00877 dir = VectorDirection(&tempX,&tempY,&tempZ); 00878 VectorOrder(dir,posinic,&tempX,&tempY,&tempZ); 00879 00880 for(j=0; j<(*Sizes)[i]; j++) 00881 { 00882 lstX.push_back(tempX[j]); 00883 lstY.push_back(tempY[j]); 00884 lstZ.push_back((*InZ)[j+cont]); 00885 } 00886 cont = cont + (*Sizes)[i]; 00887 } 00888 00889 //Fill the Finally lst in X,Y,Z --------------- 00890 int sizetemp = lstX.size(); 00891 //printf("\nJSTG-PropContour::PreparePointsForSpline"); 00892 InX->clear(); 00893 InY->clear(); 00894 InZ->clear(); 00895 for(i=0; i<sizetemp; i++) 00896 { 00897 //printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]); 00898 InX->push_back(lstX[i]); 00899 InY->push_back(lstY[i]); 00900 InZ->push_back(lstZ[i]); 00901 } 00902 00903 00904 } 00905 //---------------------------------------------------------------------------------------------------- 00906 vtkImageData* PropContour::method_Spline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ, std::vector<int>*Sizes) 00907 { 00908 //long interval = wxGetElapsedTime(TRUE); // JPRx 00909 00910 int i,j,k,sizeX,sizeS,sizeInS; 00911 int numspline; 00912 double x,y,z; 00913 //EED double spc[3]={1,1,1}; 00914 std::vector<double> lstX; 00915 std::vector<double> lstY; 00916 std::vector<double> lstZ; 00917 std::vector<double> tempX; 00918 std::vector<double> tempY; 00919 std::vector<double> tempZ; 00920 std::vector<double> tempS; 00921 _mContourModel = new manualContourModel(); 00922 _mContourModel->SetNumberOfPointsSpline(_interpnumber); 00923 00924 imagedataValue=NULL; 00925 //EED 00926 // _dimImage[0] = 200; // X axis 00927 // _dimImage[1] = 200; // Y axis 00928 // _dimImage[2] = 200; // Z axis 00929 // 00930 // unsigned short *pValue; 00931 // imagedataValue = vtkImageData::New(); 00932 // imagedataValue->SetScalarTypeToUnsignedShort(); 00933 // imagedataValue->SetSpacing(spc); 00934 // imagedataValue->SetDimensions(_dimImage); 00935 // imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1); 00936 // imagedataValue->AllocateScalars(); 00937 // imagedataValue->Update(); 00938 00939 // lstX.clear(); 00940 // lstY.clear(); 00941 // lstZ.clear(); 00942 00943 sizeX = InX->size(); 00944 for(i=0; i<sizeX; i++) 00945 { 00946 lstX.push_back((*InX)[i]); 00947 lstY.push_back((*InY)[i]); 00948 lstZ.push_back((*InZ)[i]); 00949 } 00950 00951 PreparePointsForSpline(&lstX,&lstY,&lstZ,Sizes); 00952 /*int sizetemp = lstX.size(); 00953 printf("\nJSTG-PropContour::method_Spline"); 00954 for(i=0; i<sizetemp; i++) 00955 { 00956 printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]); 00957 }*/ 00958 00959 sizeS = Sizes->size(); 00960 int cont = 0; 00961 for(i=0; i<sizeS; i++) 00962 { 00963 _mContourModel->DeleteAllPoints(); 00964 _mContourModel->SetCloseContour(true); 00965 sizeInS = (*Sizes)[i]; 00966 for(j=0; j<sizeInS; j++) 00967 { 00968 _mContourModel->AddPoint(lstX[cont],lstY[cont],lstZ[cont]); 00969 cont ++; 00970 } 00971 _mContourModel->UpdateSpline(); 00972 numspline = _mContourModel->GetNumberOfPointsSpline(); 00973 for(k=0; k<numspline; k++) 00974 { 00975 _mContourModel->GetSpline_i_Point(k,&x,&y,&z); 00976 //EED pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z); 00977 // *pValue = 128; 00978 tempX.push_back(x); 00979 tempY.push_back(y); 00980 tempZ.push_back(z); 00981 } 00982 } 00983 00984 int tam = numspline; 00985 std::vector<Vector>::iterator it; 00986 00987 ResetPlaneVector(); 00988 00989 for(i=0; i<numspline; i++) 00990 { 00991 Vector *vec = new Vector(); 00992 _planevector.push_back(*vec); 00993 }// for i 00994 00995 for(j=0; j<tam; j++) 00996 { 00997 _mContourModel->DeleteAllPoints(); 00998 _mContourModel->SetCloseContour(false); 00999 cont = 0; 01000 for(i=0; i<sizeS; i++) 01001 { 01002 //double hh=tempZ[cont+j]; // JPRx 01003 _mContourModel->AddPoint(tempX[cont+j],tempY[cont+j],tempZ[cont+j]); 01004 cont = cont + tam; 01005 } // for i 01006 _mContourModel->UpdateSpline(); 01007 numspline = _mContourModel->GetNumberOfPointsSpline(); 01008 for(k=0; k<numspline; k++) 01009 { 01010 //EED002 01011 _mContourModel->GetSpline_i_Point(k,&x,&y,&z); 01012 _planevector[k].set_x(x); 01013 _planevector[k].set_y(y); 01014 _planevector[k].set_z(z); 01015 _planevector[k].set_plane(k); 01016 //EED pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z); 01017 //EED *pValue = 128; 01018 } // for k 01019 } // for j 01020 01021 /* 01022 int tempsize = _planevector.size(); 01023 int sizexx; 01024 for(i=0; i<tempsize; i++) 01025 { 01026 sizexx = _planevector[i].getsize_x(); 01027 } 01028 */ 01029 /* 01030 interval = wxGetElapsedTime(FALSE); 01031 int zplanes = abs(tempZ[0]-tempZ[tempZ.size()-1]); 01032 long interPlane = interval/zplanes; 01033 int pointsSize = 0; 01034 for(i=0; i<sizeS; i++) 01035 { 01036 pointsSize = pointsSize + (*Sizes)[i]; 01037 } 01038 for(i=0; i<pointsSize; i++) 01039 { 01040 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*InX)[i] , (int)(*InY)[i] ,(int)(*InZ)[i]); 01041 *pValue=255; 01042 } 01043 01044 printf("\n\n JSTG - PropContour::method_Spline ------------------------"); 01045 printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %d (ms)",interval); 01046 //printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC); 01047 printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %d (ms)",interPlane); 01048 printf("\n NUMBER OF PLANES................................. %d",zplanes); 01049 printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize); 01050 printf("\n ------------------------------------------------------------"); 01051 */ 01052 return imagedataValue; 01053 } 01054 01055 //--------------------------------------------------------------------------------------------------- 01056 void PropContour::ResetPlaneVector() 01057 { 01058 // Vector *vec; 01059 // int ii,iiSize = _planevector.size(); 01060 // for (ii=0 ; ii<iiSize ; ii++) 01061 // { 01062 // vec = &(_planevector[ii]); 01063 // delete vec; 01064 // } 01065 _planevector.clear(); 01066 } 01067 01068 01069 //--------------------------------------------------------------------------------------------------- 01070 void PropContour::ResetKeyContours() 01071 { 01072 _KeyContourSizes.clear(); 01073 _KeyContourX.clear(); 01074 _KeyContourY.clear(); 01075 _KeyContourZ.clear(); 01076 } 01077 //--------------------------------------------------------------------------------------------------- 01078 void PropContour::SetKeyContours(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ) 01079 { 01080 01081 int idKeyContour = 0; 01082 int idKeyContourSizes = 0; 01083 int tmpIdKeyContSizes = 0; 01084 bool okFind = false; 01085 int i; 01086 int sizeKeyContour,Z=(int)(*InZ)[0]; 01087 sizeKeyContour = _KeyContourZ.size(); 01088 for (i=0; i<sizeKeyContour; i++) 01089 { 01090 if (i>0) 01091 { 01092 if ( (_KeyContourZ[i-1]<Z) && (_KeyContourZ[i]>=Z) ) 01093 { 01094 idKeyContour = i; 01095 idKeyContourSizes = tmpIdKeyContSizes; 01096 okFind=true; 01097 i=sizeKeyContour; 01098 } 01099 if ( (i<sizeKeyContour) && (_KeyContourZ[i-1] != _KeyContourZ[i]) ) 01100 { 01101 tmpIdKeyContSizes++; 01102 } 01103 } else { 01104 if (_KeyContourZ[0]>Z) 01105 { 01106 idKeyContour = 0; 01107 idKeyContourSizes = 0; 01108 okFind = true; 01109 i = sizeKeyContour; 01110 } 01111 } // if >O 01112 } // for 01113 01114 if (okFind==false) 01115 { 01116 idKeyContour = _KeyContourX.size(); 01117 idKeyContourSizes = _KeyContourSizes.size(); 01118 okFind=true; 01119 } 01120 01121 _KeyContourSizes.insert( _KeyContourSizes.begin() + idKeyContourSizes , InX->size() ); 01122 for(i=0; i<(int)(InX->size()); i++) 01123 { 01124 _KeyContourX.insert( _KeyContourX.begin() + idKeyContour, (*InX)[i] ); 01125 _KeyContourY.insert( _KeyContourY.begin() + idKeyContour, (*InY)[i] ); 01126 _KeyContourZ.insert( _KeyContourZ.begin() + idKeyContour, (*InZ)[i] ); 01127 } 01128 01129 01130 //EED 01131 // _KeyContourSizes.push_back( InX->size() ); 01132 // for(i=0; i<InX->size(); i++) 01133 // { 01134 // _KeyContourX.push_back( (*InX)[i] ); 01135 // _KeyContourY.push_back( (*InY)[i] ); 01136 // _KeyContourZ.push_back( (*InZ)[i] ); 01137 // } 01138 01139 } 01140 //--------------------------------------------------------------------------------------------------- 01141 vtkImageData* PropContour::CalculeSplinePropagation() 01142 { 01143 if(_KeyContourSizes.size() <= 0) 01144 { 01145 printf("\n There would be at last 1 contour"); 01146 return NULL; 01147 } 01148 if(_KeyContourSizes.size() == 1) 01149 { 01150 return NULL; 01151 } 01152 if(_KeyContourSizes.size() >= 2) 01153 { 01154 return method_Spline(&_KeyContourX,&_KeyContourY,&_KeyContourZ,&_KeyContourSizes); 01155 } 01156 return NULL; 01157 } 01158 //--------------------------------------------------------------------------------------------------- 01159 void PropContour::GetKeyContours(std::vector<double>*KeyX, std::vector<double>*KeyY, std::vector<double>*KeyZ, std::vector<int>*KeyS) 01160 { 01161 int i; 01162 KeyX->clear(); 01163 KeyY->clear(); 01164 KeyZ->clear(); 01165 KeyS->clear(); 01166 01167 for(i=0; i<(int)(_KeyContourSizes.size()); i++) 01168 { 01169 KeyS->push_back( _KeyContourSizes[i] ); 01170 } 01171 for(i=0; i<(int)(_KeyContourX.size()); i++) 01172 { 01173 KeyX->push_back( _KeyContourX[i] ); 01174 KeyY->push_back( _KeyContourY[i] ); 01175 KeyZ->push_back( _KeyContourZ[i] ); 01176 } 01177 } 01178 //--------------------------------------------------------------------------------------------------- 01179 void PropContour::GetPropagatedContours( std::vector<Vector>*planevec ) 01180 { 01181 int i,j; 01182 planevec->clear(); 01183 for(i=0; i<(int)(_planevector.size()); i++) 01184 { 01185 Vector *temp = new Vector(); 01186 temp->set_plane( _planevector[i].get_plane() ); 01187 for(j=0; j<_planevector[i].getsize_x(); j++) 01188 { 01189 temp->set_x( _planevector[i].get_x(j) ); 01190 temp->set_y( _planevector[i].get_y(j) ); 01191 temp->set_z( _planevector[i].get_z(j) ); 01192 } 01193 planevec->push_back(*temp); 01194 delete temp; 01195 } 01196 } 01197 //--------------------------------------------------------------------------------------------------- 01198 void PropContour::SetInterpNumber(int val) 01199 { 01200 _interpnumber = val; 01201 } 01202 01203 int PropContour::FindIdWithZ(double z) 01204 { 01205 int result_ID=0; 01206 int k,size=_planevector.size(); 01207 double dist,minDist=99999999; 01208 for ( k=0 ; k<size ; k++) 01209 { 01210 dist=fabs( z - _planevector[k].get_z(0) ); 01211 if (dist<minDist) 01212 { 01213 minDist=dist; 01214 result_ID = k; 01215 } 01216 }// for i 01217 return result_ID; 01218 } 01219 01220 //--------------------------------------------------------------------------------------------------- 01221 void PropContour::GetIdContour(int id, std::vector<double>*vecX, std::vector<double>*vecY, std::vector<double>*vecZ) 01222 { 01223 int i; 01224 vecX->clear(); 01225 vecY->clear(); 01226 vecZ->clear(); 01227 //int sizeplane = _planevector[id].getsize_x(); 01228 double tempx; 01229 for(i=0; i<_planevector[id].getsize_x(); i++) 01230 { 01231 vecX->push_back( _planevector[id].get_x(i) ); 01232 tempx = _planevector[id].get_x(i); 01233 vecY->push_back( _planevector[id].get_y(i) ); 01234 vecZ->push_back( _planevector[id].get_z(i) ); 01235 } 01236 } 01237 //--------------------------------------------------------------------------------------------------- 01238 //--------------------------------------------------------------------------------------------------- 01239 //---------------------------------------------------------------------------------------------------
1.6.1