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);
730 for(int i=0; i<size; i++)
732 fscanf(fd,"%lf %lf %d",&x,&y,&z);
737 SetKeyContours(&tempX,&tempY,&tempZ);
743 //---------------------------------------------------------------------------------------------------------
744 int PropContour::VectorDirection(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
747 double SumX = 0,SumY = 0;
748 double ax,ay,bx,by,axb;
749 int size = InX->size();
750 for(i=0; i<size; i++)
752 SumX = SumX + (*InX)[i];
753 SumY = SumY + (*InY)[i];
755 SumX = SumX/size; //Mass Center: X coord
756 SumY = SumY/size; //Mass Center: Y coord
760 for(i=0; i<size; i++)
764 bx = (*InX)[i+1]-SumX;
765 by = (*InY)[i+1]-SumY;
766 axb = (ax*by) - (bx*ay);
776 if(positive >= negative)
787 //----------------------------------------------------------------------------------------------------
788 void PropContour::VectorOrder(int dir, int posini, std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
790 int size = InX->size();
792 std::vector<double> tempX;
793 std::vector<double> tempY;
794 std::vector<double> tempZ;
800 for(i=0; i<size; i++)
804 tempX.push_back((*InX)[posini]);
805 tempY.push_back((*InY)[posini]);
806 tempZ.push_back((*InZ)[posini]);
815 tempX.push_back((*InX)[posini]);
816 tempY.push_back((*InY)[posini]);
817 tempZ.push_back((*InZ)[posini]);
828 for(i=0; i<size; i++)
830 InX->push_back(tempX[i]);
831 InY->push_back(tempY[i]);
832 InZ->push_back(tempZ[i]);
836 //----------------------------------------------------------------------------------------------------
837 void PropContour::PreparePointsForSpline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ,
838 std::vector<int>*Sizes)
840 int sizeS = Sizes->size();
841 int sizeV = InX->size();
842 int i,j,mem,posinic,dir,cont;
845 std::vector<double> tempX;
846 std::vector<double> tempY;
847 std::vector<double> tempZ;
848 std::vector<double> lstX;
849 std::vector<double> lstY;
850 std::vector<double> lstZ;
858 for(i=0; i<sizeS; i++)
864 for(j=0; j<(*Sizes)[i]; j++)
866 tempX.push_back((*InX)[j+mem]);
867 tempY.push_back((*InY)[j+mem]);
868 tempZ.push_back((*InZ)[j+mem]);
869 if( (*InX)[j] < leX )
875 mem = mem + (*Sizes)[i];
876 dir = VectorDirection(&tempX,&tempY,&tempZ);
877 VectorOrder(dir,posinic,&tempX,&tempY,&tempZ);
879 for(j=0; j<(*Sizes)[i]; j++)
881 lstX.push_back(tempX[j]);
882 lstY.push_back(tempY[j]);
883 lstZ.push_back((*InZ)[j+cont]);
885 cont = cont + (*Sizes)[i];
888 //Fill the Finally lst in X,Y,Z ---------------
889 int sizetemp = lstX.size();
890 //printf("\nJSTG-PropContour::PreparePointsForSpline");
894 for(i=0; i<sizetemp; i++)
896 //printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]);
897 InX->push_back(lstX[i]);
898 InY->push_back(lstY[i]);
899 InZ->push_back(lstZ[i]);
904 //----------------------------------------------------------------------------------------------------
905 vtkImageData* PropContour::method_Spline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ, std::vector<int>*Sizes)
907 long interval = wxGetElapsedTime(TRUE);
909 int i,j,k,sizeX,sizeS,sizeInS;
912 //EED double spc[3]={1,1,1};
913 std::vector<double> lstX;
914 std::vector<double> lstY;
915 std::vector<double> lstZ;
916 std::vector<double> tempX;
917 std::vector<double> tempY;
918 std::vector<double> tempZ;
919 std::vector<double> tempS;
920 _mContourModel = new manualContourModel();
921 _mContourModel->SetNumberOfPointsSpline(_interpnumber);
925 // _dimImage[0] = 200; // X axis
926 // _dimImage[1] = 200; // Y axis
927 // _dimImage[2] = 200; // Z axis
929 // unsigned short *pValue;
930 // imagedataValue = vtkImageData::New();
931 // imagedataValue->SetScalarTypeToUnsignedShort();
932 // imagedataValue->SetSpacing(spc);
933 // imagedataValue->SetDimensions(_dimImage);
934 // imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
935 // imagedataValue->AllocateScalars();
936 // imagedataValue->Update();
943 for(i=0; i<sizeX; i++)
945 lstX.push_back((*InX)[i]);
946 lstY.push_back((*InY)[i]);
947 lstZ.push_back((*InZ)[i]);
950 PreparePointsForSpline(&lstX,&lstY,&lstZ,Sizes);
951 /*int sizetemp = lstX.size();
952 printf("\nJSTG-PropContour::method_Spline");
953 for(i=0; i<sizetemp; i++)
955 printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]);
958 sizeS = Sizes->size();
960 for(i=0; i<sizeS; i++)
962 _mContourModel->DeleteAllPoints();
963 _mContourModel->SetCloseContour(true);
964 sizeInS = (*Sizes)[i];
965 for(j=0; j<sizeInS; j++)
967 _mContourModel->AddPoint(lstX[cont],lstY[cont],lstZ[cont]);
970 _mContourModel->UpdateSpline();
971 numspline = _mContourModel->GetNumberOfPointsSpline();
972 for(k=0; k<numspline; k++)
974 _mContourModel->GetSpline_i_Point(k,&x,&y,&z);
975 //EED pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z);
984 std::vector<Vector>::iterator it;
988 for(i=0; i<numspline; i++)
990 Vector *vec = new Vector();
991 _planevector.push_back(*vec);
996 _mContourModel->DeleteAllPoints();
997 _mContourModel->SetCloseContour(false);
999 for(i=0; i<sizeS; i++)
1001 double hh=tempZ[cont+j];
1002 _mContourModel->AddPoint(tempX[cont+j],tempY[cont+j],tempZ[cont+j]);
1005 _mContourModel->UpdateSpline();
1006 numspline = _mContourModel->GetNumberOfPointsSpline();
1007 for(k=0; k<numspline; k++)
1010 _mContourModel->GetSpline_i_Point(k,&x,&y,&z);
1011 _planevector[k].set_x(x);
1012 _planevector[k].set_y(y);
1013 _planevector[k].set_z(z);
1014 _planevector[k].set_plane(k);
1015 //EED pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z);
1016 //EED *pValue = 128;
1021 int tempsize = _planevector.size();
1023 for(i=0; i<tempsize; i++)
1025 sizexx = _planevector[i].getsize_x();
1029 interval = wxGetElapsedTime(FALSE);
1030 int zplanes = abs(tempZ[0]-tempZ[tempZ.size()-1]);
1031 long interPlane = interval/zplanes;
1033 for(i=0; i<sizeS; i++)
1035 pointsSize = pointsSize + (*Sizes)[i];
1037 for(i=0; i<pointsSize; i++)
1039 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*InX)[i] , (int)(*InY)[i] ,(int)(*InZ)[i]);
1043 printf("\n\n JSTG - PropContour::method_Spline ------------------------");
1044 printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %d (ms)",interval);
1045 //printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
1046 printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %d (ms)",interPlane);
1047 printf("\n NUMBER OF PLANES................................. %d",zplanes);
1048 printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
1049 printf("\n ------------------------------------------------------------");
1051 return imagedataValue;
1054 //---------------------------------------------------------------------------------------------------
1055 void PropContour::ResetPlaneVector()
1058 // int ii,iiSize = _planevector.size();
1059 // for (ii=0 ; ii<iiSize ; ii++)
1061 // vec = &(_planevector[ii]);
1064 _planevector.clear();
1068 //---------------------------------------------------------------------------------------------------
1069 void PropContour::ResetKeyContours()
1071 _KeyContourSizes.clear();
1072 _KeyContourX.clear();
1073 _KeyContourY.clear();
1074 _KeyContourZ.clear();
1076 //---------------------------------------------------------------------------------------------------
1077 void PropContour::SetKeyContours(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
1080 int idKeyContour = 0;
1081 int idKeyContourSizes = 0;
1082 int tmpIdKeyContSizes = 0;
1083 bool okFind = false;
1085 int sizeKeyContour,Z=(*InZ)[0];
1086 sizeKeyContour = _KeyContourZ.size();
1087 for (i=0; i<sizeKeyContour; i++)
1091 if ( (_KeyContourZ[i-1]<Z) && (_KeyContourZ[i]>=Z) )
1094 idKeyContourSizes = tmpIdKeyContSizes;
1098 if ( (_KeyContourZ[i-1] != _KeyContourZ[i]) )
1100 tmpIdKeyContSizes++;
1103 if (_KeyContourZ[0]>Z)
1106 idKeyContourSizes = 0;
1115 idKeyContour = _KeyContourX.size();
1116 idKeyContourSizes = _KeyContourSizes.size();
1120 _KeyContourSizes.insert( _KeyContourSizes.begin() + idKeyContourSizes , InX->size() );
1121 for(i=0; i<InX->size(); i++)
1123 _KeyContourX.insert( _KeyContourX.begin() + idKeyContour, (*InX)[i] );
1124 _KeyContourY.insert( _KeyContourY.begin() + idKeyContour, (*InY)[i] );
1125 _KeyContourZ.insert( _KeyContourZ.begin() + idKeyContour, (*InZ)[i] );
1130 // _KeyContourSizes.push_back( InX->size() );
1131 // for(i=0; i<InX->size(); i++)
1133 // _KeyContourX.push_back( (*InX)[i] );
1134 // _KeyContourY.push_back( (*InY)[i] );
1135 // _KeyContourZ.push_back( (*InZ)[i] );
1139 //---------------------------------------------------------------------------------------------------
1140 vtkImageData* PropContour::CalculeSplinePropagation()
1142 if(_KeyContourSizes.size() <= 0)
1144 printf("\n There would be at last 1 contour");
1147 if(_KeyContourSizes.size() == 1)
1151 if(_KeyContourSizes.size() >= 2)
1153 return method_Spline(&_KeyContourX,&_KeyContourY,&_KeyContourZ,&_KeyContourSizes);
1156 //---------------------------------------------------------------------------------------------------
1157 void PropContour::GetKeyContours(std::vector<double>*KeyX, std::vector<double>*KeyY, std::vector<double>*KeyZ, std::vector<int>*KeyS)
1165 for(i=0; i<_KeyContourSizes.size(); i++)
1167 KeyS->push_back( _KeyContourSizes[i] );
1169 for(i=0; i<_KeyContourX.size(); i++)
1171 KeyX->push_back( _KeyContourX[i] );
1172 KeyY->push_back( _KeyContourY[i] );
1173 KeyZ->push_back( _KeyContourZ[i] );
1176 //---------------------------------------------------------------------------------------------------
1177 void PropContour::GetPropagatedContours( std::vector<Vector>*planevec )
1181 for(i=0; i<_planevector.size(); i++)
1183 Vector *temp = new Vector();
1184 temp->set_plane( _planevector[i].get_plane() );
1185 for(j=0; j<_planevector[i].getsize_x(); j++)
1187 temp->set_x( _planevector[i].get_x(j) );
1188 temp->set_y( _planevector[i].get_y(j) );
1189 temp->set_z( _planevector[i].get_z(j) );
1191 planevec->push_back(*temp);
1195 //---------------------------------------------------------------------------------------------------
1196 void PropContour::SetInterpNumber(int val)
1198 _interpnumber = val;
1201 int PropContour::FindIdWithZ(double z)
1204 int k,size=_planevector.size();
1205 double dist,minDist=99999999;
1206 for ( k=0 ; k<size ; k++)
1208 dist=fabs( z - _planevector[k].get_z(0) );
1218 //---------------------------------------------------------------------------------------------------
1219 void PropContour::GetIdContour(int id, std::vector<double>*vecX, std::vector<double>*vecY, std::vector<double>*vecZ)
1225 int sizeplane = _planevector[id].getsize_x();
1227 for(i=0; i<_planevector[id].getsize_x(); i++)
1229 vecX->push_back( _planevector[id].get_x(i) );
1230 tempx = _planevector[id].get_x(i);
1231 vecY->push_back( _planevector[id].get_y(i) );
1232 vecZ->push_back( _planevector[id].get_z(i) );
1235 //---------------------------------------------------------------------------------------------------
1236 //---------------------------------------------------------------------------------------------------
1237 //---------------------------------------------------------------------------------------------------