1 #include "Propagation.h"
3 //------------------------------------------------------------------------------------------------------------------------------------------
4 //CLASS: Vector ------------------------------------------------------------------------------------------------------------------------
5 //------------------------------------------------------------------------------------------------------------------------------------------
21 //------------------------------------------------------------------------------------------------------------------------------------------
22 void Vector::set_vec(double val)
26 //------------------------------------------------------------------------------------------------------------------------------------------
27 void Vector::set_var(double val)
31 //------------------------------------------------------------------------------------------------------------------------------------------
32 double Vector::get_vec(int id)
40 //------------------------------------------------------------------------------------------------------------------------------------------
41 double Vector::get_var()
45 //------------------------------------------------------------------------------------------------------------------------------------------
46 int Vector::getsize_vec()
50 //------------------------------------------------------------------------------------------------------------------------------------------
51 void Vector::copyVector( std::vector<Vector>*vec1, std::vector<Vector>*vec2 )
54 if( vec1->size() != 0 )
57 for(i=0; i<vec1->size(); i++)
59 Vector *temp = new Vector();
60 temp->set_var( (*vec1)[i].get_var() );
61 for(j=0; j<(*vec1)[i].getsize_vec(); j++)
63 temp->set_vec( (*vec1)[i].get_vec(j) );
65 vec2->push_back(*temp);
70 //------------------------------------------------------------------------------------------------------------------------------------------
71 void Vector::printVector(std::vector<Vector>*vec1)
74 for(i=0; i<vec1->size(); i++)
76 printf("\n Pos (%d) => var = %f",i,(*vec1)[i].get_var());
77 for(j=0; j<(*vec1)[i].getsize_vec(); j++)
79 printf("\n vec(%d) = %f",j,(*vec1)[i].get_vec(j));
83 //------------------------------------------------------------------------------------------------------------------------------------------
84 void Vector::set_plane(int val)
88 //------------------------------------------------------------------------------------------------------------------------------------------
89 void Vector::set_x(double val)
93 //------------------------------------------------------------------------------------------------------------------------------------------
94 //------------------------------------------------------------------------------------------------------------------------------------------
95 void Vector::set_y(double val)
99 //------------------------------------------------------------------------------------------------------------------------------------------
100 //------------------------------------------------------------------------------------------------------------------------------------------
101 void Vector::set_z(double val)
103 _vecZ.push_back(val);
105 //------------------------------------------------------------------------------------------------------------------------------------------
106 int Vector::get_plane()
110 //------------------------------------------------------------------------------------------------------------------------------------------
111 double Vector::get_x(int id)
113 if( (-1<id) && (id<_vecX.size()) )
119 //------------------------------------------------------------------------------------------------------------------------------------------
120 int Vector::getsize_x()
122 if(_vecX.size() != 0)
128 //------------------------------------------------------------------------------------------------------------------------------------------
129 double Vector::get_y(int id)
131 if( (-1<id) && (id<_vecY.size()) )
137 //------------------------------------------------------------------------------------------------------------------------------------------
138 int Vector::getsize_y()
140 if(_vecY.size() != 0)
146 //------------------------------------------------------------------------------------------------------------------------------------------
147 int Vector::getsize_z()
149 if(_vecZ.size() != 0)
155 //------------------------------------------------------------------------------------------------------------------------------------------
156 double Vector::get_z(int id)
158 if( (-1<id) && (id<_vecZ.size()) )
164 //------------------------------------------------------------------------------------------------------------------------------------------
165 std::vector<double> Vector::getVec()
169 //------------------------------------------------------------------------------------------------------------------------------------------
170 void Vector::resetVec()
174 //------------------------------------------------------------------------------------------------------------------------------------------
175 //------------------------------------------------------------------------------------------------------------------------------------------
176 //------------------------------------------------------------------------------------------------------------------------------------------
180 PropContour::PropContour()
186 PropContour::~PropContour()
193 //------------------------------------------------------------------------------------
194 double PropContour::RBF_WendLand(double norm, double m_rad)
198 y = pow( 1-norm,4 ) * ( (4*norm) + 1 );
205 //------------------------------------------------------------------------------------
206 double PropContour::RBF_ThinPlate(double norm)
215 y = pow(norm,2)*log(norm);
219 //------------------------------------------------------------------------------------
220 vtkImageData* PropContour::method_RBF ( double rad, std::vector<double>*CoordX, std::vector<double>*CoordY,
221 std::vector<double>*CoordZ )
223 _dimImage[0] = 100; // X Axis
224 _dimImage[1] = 100; // Y Axis
225 _dimImage[2] = 1; // Z Axis
226 int pointsSize = CoordX->size();
227 double spc[3]={1,1,1};
228 double norm = 0, val = 0;
230 vnl_matrix<double> A( pointsSize,pointsSize, 0.0 );
231 vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 );
232 vnl_matrix<double> I( pointsSize,pointsSize, 0.0 );
233 vnl_vector<double> H( pointsSize, 1 );
234 vnl_vector<double> D( pointsSize, 0.0 );
235 vnl_vector<double> Q( 2, 0.0 );
236 vnl_matrix<double> Impl( _dimImage[0],_dimImage[1], 0.0 );
238 unsigned short *pValue;
239 imagedataValue = vtkImageData::New();
240 imagedataValue->SetScalarTypeToUnsignedShort();
241 imagedataValue->SetSpacing(spc);
242 imagedataValue->SetDimensions(_dimImage);
243 imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
244 imagedataValue->AllocateScalars();
245 imagedataValue->Update();
248 for(i=0; i<pointsSize; i++)
250 for(j=0; j<pointsSize; j++)
252 norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) );
253 A(i,j) = RBF_WendLand(norm,rad);
254 //A(i,j) = RBF_ThinPlate(norm);
258 Ainv = vnl_matrix_inverse<double>(A);
261 for(i=0; i < _dimImage[1]; i++)
263 for(j=0; j < _dimImage[0]; j++)
268 for(h=0; h<pointsSize; h++)
270 norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) );
271 //val = val + ( D(h) * RBF_ThinPlate(norm) );
272 val = val + ( D(h) * RBF_WendLand(norm,rad) );
274 Impl(i,j) = val - 1.0;
276 // if ( (Impl(i,j)>=-0.001) && (Impl(i,j) <=0.001) )
283 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,0);
284 // vvalue = (Impl(i,j)*256+256);
292 pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,0);
304 for(i=0; i<pointsSize; i=i+5)
306 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,0);
310 return imagedataValue;
312 //----------------------------------------------------------------------------------------------------
313 double PropContour::RBF_ThinPlate_3D(double norm)
315 return norm = pow( norm,3 );
317 //----------------------------------------------------------------------------------------------------
318 vtkImageData* PropContour::method_RBF_3D ( double rad, std::vector<double>*CoordX, std::vector<double>*CoordY,
319 std::vector<double>*CoordZ )
321 long interval = wxGetElapsedTime(TRUE);
324 int pointsSize = CoordX->size();
325 double max = -1, min = 1000, minz = 100000, maxz = -10000;
327 for( i=0; i<pointsSize; i++ )
329 if( (*CoordX)[i] > max )
333 if( (*CoordY)[i] > max )
337 if( (*CoordX)[i] < min )
341 if( (*CoordY)[i] < min )
345 if( (*CoordZ)[i] < minz )
349 if( (*CoordZ)[i] > maxz )
355 _dimImage[0] = 200; // X axis
356 _dimImage[1] = 200; // Y axis
357 _dimImage[2] = 200; // Z axis
359 double spc[3]={1,1,1};
360 double norm = 0, val, vvalue;
362 vnl_matrix<double> A( pointsSize,pointsSize, 0.0 );
363 vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 );
364 vnl_matrix<double> I( pointsSize,pointsSize, 0.0 );
365 vnl_vector<double> H( pointsSize, 1 );
366 vnl_vector<double> D( pointsSize, 0.0 );
367 vnl_vector<double> Q( 3, 0.0 );
369 unsigned short *pValue;
370 imagedataValue = vtkImageData::New();
371 imagedataValue->SetScalarTypeToUnsignedShort();
372 imagedataValue->SetSpacing(spc);
373 imagedataValue->SetDimensions(_dimImage);
374 imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
375 imagedataValue->AllocateScalars();
376 imagedataValue->Update();
378 for(i=0; i<pointsSize; i++)
380 for(j=0; j<pointsSize; j++)
382 norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) + pow((*CoordZ)[i]-(*CoordZ)[j],2) );
383 A(i,j) = RBF_WendLand(norm,rad);
386 D = vnl_matrix_inverse<double>(A)* H;
389 // for(k=0; k<_dimImage[2]; k++)
391 // for(j=0; j<_dimImage[1]; j++)
393 // for(i=0; i<_dimImage[0]; i++)
395 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
402 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(0,0,0);
406 // for(h=0; h<pointsSize; h++)
408 // norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
409 // //val = val + ( D(h) * RBF_ThinPlate_3D(norm) );
410 // val = val + ( D(h) * RBF_WendLand(norm,rad) );
411 // if( h == pointsSize-1 )
414 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
427 // if( j == _dimImage[0] )
432 // if( i == _dimImage[1])
438 // if(k == _dimImage[2])
452 for(h=0; h<pointsSize; h++)
454 norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
455 val = val + ( D(h) * RBF_WendLand(norm,rad) );
456 if( h == pointsSize-1 )
459 pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
472 if( j == (int)max+10)
477 if( i == (int)max+10)
492 interval = wxGetElapsedTime(FALSE);
493 long interPlane = interval/_dimImage[2];
495 for(i=0; i<pointsSize; i++)
497 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,(int)(*CoordZ)[i]);
501 long intervalPC = wxGetElapsedTime();
503 printf("\n\n JSTG - PropContour::method_RBF_3D ------------------------");
504 printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %lld (ms)",interval);
505 printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
506 printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %lld (ms)",interPlane);
507 printf("\n NUMBER OF PLANES................................. %d",k);
508 printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
509 printf("\n ------------------------------------------------------------");
511 return imagedataValue;
513 //----------------------------------------------------------------------------------------------------
515 vtkImageData* PropContour::method_RBF_3D_ThinPlate ( double rad, std::vector<double>*CoordX, std::vector<double>*CoordY,
516 std::vector<double>*CoordZ )
518 long interval = wxGetElapsedTime(TRUE);
521 int pointsSize = CoordX->size();
522 double max = -1, min = 1000, minz = 100000, maxz = -10000;
524 for( i=0; i<pointsSize; i++ )
526 if( (*CoordX)[i] > max )
530 if( (*CoordY)[i] > max )
534 if( (*CoordX)[i] < min )
538 if( (*CoordY)[i] < min )
542 if( (*CoordZ)[i] < minz )
546 if( (*CoordZ)[i] > maxz )
552 _dimImage[0] = 190; // X axis
553 _dimImage[1] = 190; // Y axis
554 _dimImage[2] = 190; // Z axis
556 double spc[3]={1,1,1};
557 double norm = 0, val, vvalue;
559 vnl_matrix<double> A( pointsSize,pointsSize, 0.0 );
560 vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 );
561 vnl_matrix<double> I( pointsSize,pointsSize, 0.0 );
562 vnl_vector<double> H( pointsSize, 1 );
563 vnl_vector<double> D( pointsSize, 0.0 );
564 vnl_vector<double> Q( 3, 0.0 );
566 unsigned short *pValue;
567 imagedataValue = vtkImageData::New();
568 imagedataValue->SetScalarTypeToUnsignedShort();
569 imagedataValue->SetSpacing(spc);
570 imagedataValue->SetDimensions(_dimImage);
571 imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,-10,_dimImage[2]-1);
572 imagedataValue->AllocateScalars();
573 imagedataValue->Update();
575 for(i=0; i<pointsSize; i++)
577 for(j=0; j<pointsSize; j++)
579 norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) + pow((*CoordZ)[i]-(*CoordZ)[j],2) );
580 A(i,j) = RBF_ThinPlate_3D(norm);
583 D = vnl_matrix_inverse<double>(A)* H;
586 // for(k=0; k<_dimImage[2]; k++)
588 // for(j=0; j<_dimImage[1]; j++)
590 // for(i=0; i<_dimImage[0]; i++)
592 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
600 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(0,0,0);
604 // for(h=0; h<pointsSize; h++)
606 // norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
607 // //val = val + ( D(h) * RBF_ThinPlate_3D(norm) );
608 // val = val + ( D(h) * RBF_WendLand(norm,rad) );
609 // if( h == pointsSize-1 )
612 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
625 // if( j == _dimImage[0] )
630 // if( i == _dimImage[1])
636 // if(k == _dimImage[2])
650 for(h=0; h<pointsSize; h++)
652 norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
653 val = val + ( D(h) * RBF_ThinPlate_3D(norm) );
654 if( h == pointsSize-1 )
657 pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
670 if( j == (int)max+10)
675 if( i == (int)max+10)
690 interval = wxGetElapsedTime(FALSE);
691 long interPlane = interval/_dimImage[2];
693 for(i=0; i<pointsSize; i++)
695 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,(int)(*CoordZ)[i]);
699 long intervalPC = wxGetElapsedTime();
701 printf("\n\n JSTG - PropContour::method_RBF_3D ------------------------");
702 printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %lld (ms)",interval);
703 printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
704 printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %lld (ms)",interPlane);
705 printf("\n NUMBER OF PLANES................................. %d",k);
706 printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
707 printf("\n ------------------------------------------------------------");
709 return imagedataValue;
714 //---------------------------------------------------------------------------------------------------------
715 void PropContour::ReadKeyContour(FILE* fd)
721 std::vector<double> tempX;
722 std::vector<double> tempY;
723 std::vector<double> tempZ;
729 //fscanf(fd," %s %d",&firstline,&size); // JPRx
730 fscanf(fd," %s %d",firstline,&size);
731 for(int i=0; i<size; i++)
733 fscanf(fd,"%lf %lf %d",&x,&y,&z);
738 SetKeyContours(&tempX,&tempY,&tempZ);
744 //---------------------------------------------------------------------------------------------------------
745 int PropContour::VectorDirection(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
748 double SumX = 0,SumY = 0;
749 double ax,ay,bx,by,axb;
750 int size = InX->size();
751 for(i=0; i<size; i++)
753 SumX = SumX + (*InX)[i];
754 SumY = SumY + (*InY)[i];
756 SumX = SumX/size; //Mass Center: X coord
757 SumY = SumY/size; //Mass Center: Y coord
761 for(i=0; i<size; i++)
765 bx = (*InX)[i+1]-SumX;
766 by = (*InY)[i+1]-SumY;
767 axb = (ax*by) - (bx*ay);
777 if(positive >= negative)
788 //----------------------------------------------------------------------------------------------------
789 void PropContour::VectorOrder(int dir, int posini, std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
791 int size = InX->size();
793 std::vector<double> tempX;
794 std::vector<double> tempY;
795 std::vector<double> tempZ;
801 for(i=0; i<size; i++)
805 tempX.push_back((*InX)[posini]);
806 tempY.push_back((*InY)[posini]);
807 tempZ.push_back((*InZ)[posini]);
816 tempX.push_back((*InX)[posini]);
817 tempY.push_back((*InY)[posini]);
818 tempZ.push_back((*InZ)[posini]);
829 for(i=0; i<size; i++)
831 InX->push_back(tempX[i]);
832 InY->push_back(tempY[i]);
833 InZ->push_back(tempZ[i]);
837 //----------------------------------------------------------------------------------------------------
838 void PropContour::PreparePointsForSpline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ,
839 std::vector<int>*Sizes)
841 int sizeS = Sizes->size();
842 //int sizeV = InX->size(); // JPRx
843 int i,j,mem,posinic,dir,cont;
846 std::vector<double> tempX;
847 std::vector<double> tempY;
848 std::vector<double> tempZ;
849 std::vector<double> lstX;
850 std::vector<double> lstY;
851 std::vector<double> lstZ;
859 for(i=0; i<sizeS; i++)
865 for(j=0; j<(*Sizes)[i]; j++)
867 tempX.push_back((*InX)[j+mem]);
868 tempY.push_back((*InY)[j+mem]);
869 tempZ.push_back((*InZ)[j+mem]);
870 if( (*InX)[j] < leX )
876 mem = mem + (*Sizes)[i];
877 dir = VectorDirection(&tempX,&tempY,&tempZ);
878 VectorOrder(dir,posinic,&tempX,&tempY,&tempZ);
880 for(j=0; j<(*Sizes)[i]; j++)
882 lstX.push_back(tempX[j]);
883 lstY.push_back(tempY[j]);
884 lstZ.push_back((*InZ)[j+cont]);
886 cont = cont + (*Sizes)[i];
889 //Fill the Finally lst in X,Y,Z ---------------
890 int sizetemp = lstX.size();
891 //printf("\nJSTG-PropContour::PreparePointsForSpline");
895 for(i=0; i<sizetemp; i++)
897 //printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]);
898 InX->push_back(lstX[i]);
899 InY->push_back(lstY[i]);
900 InZ->push_back(lstZ[i]);
905 //----------------------------------------------------------------------------------------------------
906 vtkImageData* PropContour::method_Spline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ, std::vector<int>*Sizes)
908 //long interval = wxGetElapsedTime(TRUE); // JPRx
910 int i,j,k,sizeX,sizeS,sizeInS;
913 //EED double spc[3]={1,1,1};
914 std::vector<double> lstX;
915 std::vector<double> lstY;
916 std::vector<double> lstZ;
917 std::vector<double> tempX;
918 std::vector<double> tempY;
919 std::vector<double> tempZ;
920 std::vector<double> tempS;
921 _mContourModel = new manualContourModel();
922 _mContourModel->SetNumberOfPointsSpline(_interpnumber);
926 // _dimImage[0] = 200; // X axis
927 // _dimImage[1] = 200; // Y axis
928 // _dimImage[2] = 200; // Z axis
930 // unsigned short *pValue;
931 // imagedataValue = vtkImageData::New();
932 // imagedataValue->SetScalarTypeToUnsignedShort();
933 // imagedataValue->SetSpacing(spc);
934 // imagedataValue->SetDimensions(_dimImage);
935 // imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
936 // imagedataValue->AllocateScalars();
937 // imagedataValue->Update();
944 for(i=0; i<sizeX; i++)
946 lstX.push_back((*InX)[i]);
947 lstY.push_back((*InY)[i]);
948 lstZ.push_back((*InZ)[i]);
951 PreparePointsForSpline(&lstX,&lstY,&lstZ,Sizes);
952 /*int sizetemp = lstX.size();
953 printf("\nJSTG-PropContour::method_Spline");
954 for(i=0; i<sizetemp; i++)
956 printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]);
959 sizeS = Sizes->size();
961 for(i=0; i<sizeS; i++)
963 _mContourModel->DeleteAllPoints();
964 _mContourModel->SetCloseContour(true);
965 sizeInS = (*Sizes)[i];
966 for(j=0; j<sizeInS; j++)
968 _mContourModel->AddPoint(lstX[cont],lstY[cont],lstZ[cont]);
971 _mContourModel->UpdateSpline();
972 numspline = _mContourModel->GetNumberOfPointsSpline();
973 for(k=0; k<numspline; k++)
975 _mContourModel->GetSpline_i_Point(k,&x,&y,&z);
976 //EED pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z);
985 std::vector<Vector>::iterator it;
989 for(i=0; i<numspline; i++)
991 Vector *vec = new Vector();
992 _planevector.push_back(*vec);
997 _mContourModel->DeleteAllPoints();
998 _mContourModel->SetCloseContour(false);
1000 for(i=0; i<sizeS; i++)
1002 //double hh=tempZ[cont+j]; // JPRx
1003 _mContourModel->AddPoint(tempX[cont+j],tempY[cont+j],tempZ[cont+j]);
1006 _mContourModel->UpdateSpline();
1007 numspline = _mContourModel->GetNumberOfPointsSpline();
1008 for(k=0; k<numspline; k++)
1011 _mContourModel->GetSpline_i_Point(k,&x,&y,&z);
1012 _planevector[k].set_x(x);
1013 _planevector[k].set_y(y);
1014 _planevector[k].set_z(z);
1015 _planevector[k].set_plane(k);
1016 //EED pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z);
1017 //EED *pValue = 128;
1022 int tempsize = _planevector.size();
1024 for(i=0; i<tempsize; i++)
1026 sizexx = _planevector[i].getsize_x();
1030 interval = wxGetElapsedTime(FALSE);
1031 int zplanes = abs(tempZ[0]-tempZ[tempZ.size()-1]);
1032 long interPlane = interval/zplanes;
1034 for(i=0; i<sizeS; i++)
1036 pointsSize = pointsSize + (*Sizes)[i];
1038 for(i=0; i<pointsSize; i++)
1040 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*InX)[i] , (int)(*InY)[i] ,(int)(*InZ)[i]);
1044 printf("\n\n JSTG - PropContour::method_Spline ------------------------");
1045 printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %d (ms)",interval);
1046 //printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
1047 printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %d (ms)",interPlane);
1048 printf("\n NUMBER OF PLANES................................. %d",zplanes);
1049 printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
1050 printf("\n ------------------------------------------------------------");
1052 return imagedataValue;
1055 //---------------------------------------------------------------------------------------------------
1056 void PropContour::ResetPlaneVector()
1059 // int ii,iiSize = _planevector.size();
1060 // for (ii=0 ; ii<iiSize ; ii++)
1062 // vec = &(_planevector[ii]);
1065 _planevector.clear();
1069 //---------------------------------------------------------------------------------------------------
1070 void PropContour::ResetKeyContours()
1072 _KeyContourSizes.clear();
1073 _KeyContourX.clear();
1074 _KeyContourY.clear();
1075 _KeyContourZ.clear();
1077 //---------------------------------------------------------------------------------------------------
1078 void PropContour::SetKeyContours(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
1081 int idKeyContour = 0;
1082 int idKeyContourSizes = 0;
1083 int tmpIdKeyContSizes = 0;
1084 bool okFind = false;
1086 int sizeKeyContour,Z=(int)(*InZ)[0];
1087 sizeKeyContour = _KeyContourZ.size();
1088 for (i=0; i<sizeKeyContour; i++)
1092 if ( (_KeyContourZ[i-1]<Z) && (_KeyContourZ[i]>=Z) )
1095 idKeyContourSizes = tmpIdKeyContSizes;
1099 if ( (_KeyContourZ[i-1] != _KeyContourZ[i]) )
1101 tmpIdKeyContSizes++;
1104 if (_KeyContourZ[0]>Z)
1107 idKeyContourSizes = 0;
1116 idKeyContour = _KeyContourX.size();
1117 idKeyContourSizes = _KeyContourSizes.size();
1121 _KeyContourSizes.insert( _KeyContourSizes.begin() + idKeyContourSizes , InX->size() );
1122 for(i=0; i<InX->size(); i++)
1124 _KeyContourX.insert( _KeyContourX.begin() + idKeyContour, (*InX)[i] );
1125 _KeyContourY.insert( _KeyContourY.begin() + idKeyContour, (*InY)[i] );
1126 _KeyContourZ.insert( _KeyContourZ.begin() + idKeyContour, (*InZ)[i] );
1131 // _KeyContourSizes.push_back( InX->size() );
1132 // for(i=0; i<InX->size(); i++)
1134 // _KeyContourX.push_back( (*InX)[i] );
1135 // _KeyContourY.push_back( (*InY)[i] );
1136 // _KeyContourZ.push_back( (*InZ)[i] );
1140 //---------------------------------------------------------------------------------------------------
1141 vtkImageData* PropContour::CalculeSplinePropagation()
1143 if(_KeyContourSizes.size() <= 0)
1145 printf("\n There would be at last 1 contour");
1148 if(_KeyContourSizes.size() == 1)
1152 if(_KeyContourSizes.size() >= 2)
1154 return method_Spline(&_KeyContourX,&_KeyContourY,&_KeyContourZ,&_KeyContourSizes);
1157 //---------------------------------------------------------------------------------------------------
1158 void PropContour::GetKeyContours(std::vector<double>*KeyX, std::vector<double>*KeyY, std::vector<double>*KeyZ, std::vector<int>*KeyS)
1166 for(i=0; i<_KeyContourSizes.size(); i++)
1168 KeyS->push_back( _KeyContourSizes[i] );
1170 for(i=0; i<_KeyContourX.size(); i++)
1172 KeyX->push_back( _KeyContourX[i] );
1173 KeyY->push_back( _KeyContourY[i] );
1174 KeyZ->push_back( _KeyContourZ[i] );
1177 //---------------------------------------------------------------------------------------------------
1178 void PropContour::GetPropagatedContours( std::vector<Vector>*planevec )
1182 for(i=0; i<_planevector.size(); i++)
1184 Vector *temp = new Vector();
1185 temp->set_plane( _planevector[i].get_plane() );
1186 for(j=0; j<_planevector[i].getsize_x(); j++)
1188 temp->set_x( _planevector[i].get_x(j) );
1189 temp->set_y( _planevector[i].get_y(j) );
1190 temp->set_z( _planevector[i].get_z(j) );
1192 planevec->push_back(*temp);
1196 //---------------------------------------------------------------------------------------------------
1197 void PropContour::SetInterpNumber(int val)
1199 _interpnumber = val;
1202 int PropContour::FindIdWithZ(double z)
1205 int k,size=_planevector.size();
1206 double dist,minDist=99999999;
1207 for ( k=0 ; k<size ; k++)
1209 dist=fabs( z - _planevector[k].get_z(0) );
1219 //---------------------------------------------------------------------------------------------------
1220 void PropContour::GetIdContour(int id, std::vector<double>*vecX, std::vector<double>*vecY, std::vector<double>*vecZ)
1226 //int sizeplane = _planevector[id].getsize_x();
1228 for(i=0; i<_planevector[id].getsize_x(); i++)
1230 vecX->push_back( _planevector[id].get_x(i) );
1231 tempx = _planevector[id].get_x(i);
1232 vecY->push_back( _planevector[id].get_y(i) );
1233 vecZ->push_back( _planevector[id].get_z(i) );
1236 //---------------------------------------------------------------------------------------------------
1237 //---------------------------------------------------------------------------------------------------
1238 //---------------------------------------------------------------------------------------------------