1 /*# ---------------------------------------------------------------------
3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
5 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
6 # Previous Authors : Laurent Guigues, Jean-Pierre Roux
7 # CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
9 # This software is governed by the CeCILL-B license under French law and
10 # abiding by the rules of distribution of free software. You can use,
11 # modify and/ or redistribute the software under the terms of the CeCILL-B
12 # license as circulated by CEA, CNRS and INRIA at the following URL
13 # http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
14 # or in the file LICENSE.txt.
16 # As a counterpart to the access to the source code and rights to copy,
17 # modify and redistribute granted by the license, users are provided only
18 # with a limited warranty and the software's author, the holder of the
19 # economic rights, and the successive licensors have only limited
22 # The fact that you are presently reading this means that you have had
23 # knowledge of the CeCILL-B license and that you accept its terms.
24 # ------------------------------------------------------------------------ */
26 #include "Propagation.h"
28 //------------------------------------------------------------------------------------------------------------------------------------------
29 //CLASS: Vector ------------------------------------------------------------------------------------------------------------------------
30 //------------------------------------------------------------------------------------------------------------------------------------------
46 //------------------------------------------------------------------------------------------------------------------------------------------
47 void Vector::set_vec(double val)
51 //------------------------------------------------------------------------------------------------------------------------------------------
52 void Vector::set_var(double val)
56 //------------------------------------------------------------------------------------------------------------------------------------------
57 double Vector::get_vec(int id)
65 //------------------------------------------------------------------------------------------------------------------------------------------
66 double Vector::get_var()
70 //------------------------------------------------------------------------------------------------------------------------------------------
71 int Vector::getsize_vec()
75 //------------------------------------------------------------------------------------------------------------------------------------------
76 void Vector::copyVector( std::vector<Vector>*vec1, std::vector<Vector>*vec2 )
79 if( vec1->size() != 0 )
82 for(i=0; i<(int)(vec1->size()); i++)
84 Vector *temp = new Vector();
85 temp->set_var( (*vec1)[i].get_var() );
86 for(j=0; j<(*vec1)[i].getsize_vec(); j++)
88 temp->set_vec( (*vec1)[i].get_vec(j) );
90 vec2->push_back(*temp);
95 //------------------------------------------------------------------------------------------------------------------------------------------
96 void Vector::printVector(std::vector<Vector>*vec1)
99 for(i=0; i<(int)(vec1->size()); i++)
101 printf("\n Pos (%d) => var = %f",i,(*vec1)[i].get_var());
102 for(j=0; j<(*vec1)[i].getsize_vec(); j++)
104 printf("\n vec(%d) = %f",j,(*vec1)[i].get_vec(j));
108 //------------------------------------------------------------------------------------------------------------------------------------------
109 void Vector::set_plane(int val)
113 //------------------------------------------------------------------------------------------------------------------------------------------
114 void Vector::set_x(double val)
116 _vecX.push_back(val);
118 //------------------------------------------------------------------------------------------------------------------------------------------
119 //------------------------------------------------------------------------------------------------------------------------------------------
120 void Vector::set_y(double val)
122 _vecY.push_back(val);
124 //------------------------------------------------------------------------------------------------------------------------------------------
125 //------------------------------------------------------------------------------------------------------------------------------------------
126 void Vector::set_z(double val)
128 _vecZ.push_back(val);
130 //------------------------------------------------------------------------------------------------------------------------------------------
131 int Vector::get_plane()
135 //------------------------------------------------------------------------------------------------------------------------------------------
136 double Vector::get_x(int id)
138 if( (-1<id) && (id<(int)(_vecX.size()) ) )
144 //------------------------------------------------------------------------------------------------------------------------------------------
145 int Vector::getsize_x()
147 if(_vecX.size() != 0)
153 //------------------------------------------------------------------------------------------------------------------------------------------
154 double Vector::get_y(int id)
156 if( (-1<id) && (id<(int)(_vecY.size())) )
162 //------------------------------------------------------------------------------------------------------------------------------------------
163 int Vector::getsize_y()
165 if(_vecY.size() != 0)
171 //------------------------------------------------------------------------------------------------------------------------------------------
172 int Vector::getsize_z()
174 if(_vecZ.size() != 0)
180 //------------------------------------------------------------------------------------------------------------------------------------------
181 double Vector::get_z(int id)
183 if( (-1<id) && (id<(int)(_vecZ.size())) )
189 //------------------------------------------------------------------------------------------------------------------------------------------
190 std::vector<double> Vector::getVec()
194 //------------------------------------------------------------------------------------------------------------------------------------------
195 void Vector::resetVec()
199 //------------------------------------------------------------------------------------------------------------------------------------------
200 //------------------------------------------------------------------------------------------------------------------------------------------
201 //------------------------------------------------------------------------------------------------------------------------------------------
205 PropContour::PropContour()
211 PropContour::~PropContour()
218 //------------------------------------------------------------------------------------
219 double PropContour::RBF_WendLand(double norm, double m_rad)
223 y = pow( 1-norm,4 ) * ( (4*norm) + 1 );
230 //------------------------------------------------------------------------------------
231 double PropContour::RBF_ThinPlate(double norm)
240 y = pow(norm,2)*log(norm);
244 //------------------------------------------------------------------------------------
245 vtkImageData* PropContour::method_RBF ( double rad, std::vector<double>*CoordX, std::vector<double>*CoordY,
246 std::vector<double>*CoordZ )
248 _dimImage[0] = 100; // X Axis
249 _dimImage[1] = 100; // Y Axis
250 _dimImage[2] = 1; // Z Axis
251 int pointsSize = CoordX->size();
252 double spc[3]={1,1,1};
253 double norm = 0, val = 0;
255 vnl_matrix<double> A( pointsSize,pointsSize, 0.0 );
256 vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 );
257 vnl_matrix<double> I( pointsSize,pointsSize, 0.0 );
258 vnl_vector<double> H( pointsSize, 1 );
259 vnl_vector<double> D( pointsSize, 0.0 );
260 vnl_vector<double> Q( 2, 0.0 );
261 vnl_matrix<double> Impl( _dimImage[0],_dimImage[1], 0.0 );
263 unsigned short *pValue;
264 imagedataValue = vtkImageData::New();
265 imagedataValue->SetScalarTypeToUnsignedShort();
266 imagedataValue->SetSpacing(spc);
267 imagedataValue->SetDimensions(_dimImage);
268 imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
269 imagedataValue->AllocateScalars();
270 imagedataValue->Update();
273 for(i=0; i<pointsSize; i++)
275 for(j=0; j<pointsSize; j++)
277 norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) );
278 A(i,j) = RBF_WendLand(norm,rad);
279 //A(i,j) = RBF_ThinPlate(norm);
283 Ainv = vnl_matrix_inverse<double>(A);
286 for(i=0; i < _dimImage[1]; i++)
288 for(j=0; j < _dimImage[0]; j++)
293 for(h=0; h<pointsSize; h++)
295 norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) );
296 //val = val + ( D(h) * RBF_ThinPlate(norm) );
297 val = val + ( D(h) * RBF_WendLand(norm,rad) );
299 Impl(i,j) = val - 1.0;
301 // if ( (Impl(i,j)>=-0.001) && (Impl(i,j) <=0.001) )
308 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,0);
309 // vvalue = (Impl(i,j)*256+256);
317 pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,0);
329 for(i=0; i<pointsSize; i=i+5)
331 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,0);
335 return imagedataValue;
337 //----------------------------------------------------------------------------------------------------
338 double PropContour::RBF_ThinPlate_3D(double norm)
340 return norm = pow( norm,3 );
342 //----------------------------------------------------------------------------------------------------
343 vtkImageData* PropContour::method_RBF_3D ( double rad, std::vector<double>*CoordX, std::vector<double>*CoordY,
344 std::vector<double>*CoordZ )
346 long interval = wxGetElapsedTime(TRUE);
349 int pointsSize = CoordX->size();
350 double max = -1, min = 1000, minz = 100000, maxz = -10000;
352 for( i=0; i<pointsSize; i++ )
354 if( (*CoordX)[i] > max )
358 if( (*CoordY)[i] > max )
362 if( (*CoordX)[i] < min )
366 if( (*CoordY)[i] < min )
370 if( (*CoordZ)[i] < minz )
374 if( (*CoordZ)[i] > maxz )
380 _dimImage[0] = 200; // X axis
381 _dimImage[1] = 200; // Y axis
382 _dimImage[2] = 200; // Z axis
384 double spc[3]={1,1,1};
385 double norm = 0, val, vvalue;
387 vnl_matrix<double> A( pointsSize,pointsSize, 0.0 );
388 vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 );
389 vnl_matrix<double> I( pointsSize,pointsSize, 0.0 );
390 vnl_vector<double> H( pointsSize, 1 );
391 vnl_vector<double> D( pointsSize, 0.0 );
392 vnl_vector<double> Q( 3, 0.0 );
394 unsigned short *pValue;
395 imagedataValue = vtkImageData::New();
396 imagedataValue->SetScalarTypeToUnsignedShort();
397 imagedataValue->SetSpacing(spc);
398 imagedataValue->SetDimensions(_dimImage);
399 imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
400 imagedataValue->AllocateScalars();
401 imagedataValue->Update();
403 for(i=0; i<pointsSize; i++)
405 for(j=0; j<pointsSize; j++)
407 norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) + pow((*CoordZ)[i]-(*CoordZ)[j],2) );
408 A(i,j) = RBF_WendLand(norm,rad);
411 D = vnl_matrix_inverse<double>(A)* H;
414 // for(k=0; k<_dimImage[2]; k++)
416 // for(j=0; j<_dimImage[1]; j++)
418 // for(i=0; i<_dimImage[0]; i++)
420 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
427 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(0,0,0);
431 // for(h=0; h<pointsSize; h++)
433 // norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
434 // //val = val + ( D(h) * RBF_ThinPlate_3D(norm) );
435 // val = val + ( D(h) * RBF_WendLand(norm,rad) );
436 // if( h == pointsSize-1 )
439 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
452 // if( j == _dimImage[0] )
457 // if( i == _dimImage[1])
463 // if(k == _dimImage[2])
477 for(h=0; h<pointsSize; h++)
479 norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
480 val = val + ( D(h) * RBF_WendLand(norm,rad) );
481 if( h == pointsSize-1 )
484 pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
497 if( j == (int)max+10)
502 if( i == (int)max+10)
517 interval = wxGetElapsedTime(FALSE);
518 long interPlane = interval/_dimImage[2];
520 for(i=0; i<pointsSize; i++)
522 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,(int)(*CoordZ)[i]);
526 long intervalPC = wxGetElapsedTime();
528 printf("\n\n JSTG - PropContour::method_RBF_3D ------------------------");
529 printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %lld (ms)",interval);
530 printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
531 printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %lld (ms)",interPlane);
532 printf("\n NUMBER OF PLANES................................. %d",k);
533 printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
534 printf("\n ------------------------------------------------------------");
536 return imagedataValue;
538 //----------------------------------------------------------------------------------------------------
540 vtkImageData* PropContour::method_RBF_3D_ThinPlate ( double rad, std::vector<double>*CoordX, std::vector<double>*CoordY,
541 std::vector<double>*CoordZ )
543 long interval = wxGetElapsedTime(TRUE);
546 int pointsSize = CoordX->size();
547 double max = -1, min = 1000, minz = 100000, maxz = -10000;
549 for( i=0; i<pointsSize; i++ )
551 if( (*CoordX)[i] > max )
555 if( (*CoordY)[i] > max )
559 if( (*CoordX)[i] < min )
563 if( (*CoordY)[i] < min )
567 if( (*CoordZ)[i] < minz )
571 if( (*CoordZ)[i] > maxz )
577 _dimImage[0] = 190; // X axis
578 _dimImage[1] = 190; // Y axis
579 _dimImage[2] = 190; // Z axis
581 double spc[3]={1,1,1};
582 double norm = 0, val, vvalue;
584 vnl_matrix<double> A( pointsSize,pointsSize, 0.0 );
585 vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 );
586 vnl_matrix<double> I( pointsSize,pointsSize, 0.0 );
587 vnl_vector<double> H( pointsSize, 1 );
588 vnl_vector<double> D( pointsSize, 0.0 );
589 vnl_vector<double> Q( 3, 0.0 );
591 unsigned short *pValue;
592 imagedataValue = vtkImageData::New();
593 imagedataValue->SetScalarTypeToUnsignedShort();
594 imagedataValue->SetSpacing(spc);
595 imagedataValue->SetDimensions(_dimImage);
596 imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,-10,_dimImage[2]-1);
597 imagedataValue->AllocateScalars();
598 imagedataValue->Update();
600 for(i=0; i<pointsSize; i++)
602 for(j=0; j<pointsSize; j++)
604 norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) + pow((*CoordZ)[i]-(*CoordZ)[j],2) );
605 A(i,j) = RBF_ThinPlate_3D(norm);
608 D = vnl_matrix_inverse<double>(A)* H;
611 // for(k=0; k<_dimImage[2]; k++)
613 // for(j=0; j<_dimImage[1]; j++)
615 // for(i=0; i<_dimImage[0]; i++)
617 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
625 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(0,0,0);
629 // for(h=0; h<pointsSize; h++)
631 // norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
632 // //val = val + ( D(h) * RBF_ThinPlate_3D(norm) );
633 // val = val + ( D(h) * RBF_WendLand(norm,rad) );
634 // if( h == pointsSize-1 )
637 // pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
650 // if( j == _dimImage[0] )
655 // if( i == _dimImage[1])
661 // if(k == _dimImage[2])
675 for(h=0; h<pointsSize; h++)
677 norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
678 val = val + ( D(h) * RBF_ThinPlate_3D(norm) );
679 if( h == pointsSize-1 )
682 pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
695 if( j == (int)max+10)
700 if( i == (int)max+10)
715 interval = wxGetElapsedTime(FALSE);
716 long interPlane = interval/_dimImage[2];
718 for(i=0; i<pointsSize; i++)
720 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,(int)(*CoordZ)[i]);
724 long intervalPC = wxGetElapsedTime();
726 printf("\n\n JSTG - PropContour::method_RBF_3D ------------------------");
727 printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %lld (ms)",interval);
728 printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
729 printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %lld (ms)",interPlane);
730 printf("\n NUMBER OF PLANES................................. %d",k);
731 printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
732 printf("\n ------------------------------------------------------------");
734 return imagedataValue;
739 //---------------------------------------------------------------------------------------------------------
740 void PropContour::ReadKeyContour(FILE* fd)
746 std::vector<double> tempX;
747 std::vector<double> tempY;
748 std::vector<double> tempZ;
754 //fscanf(fd," %s %d",&firstline,&size); // JPRx
755 fscanf(fd," %s %d",firstline,&size);
756 for(int i=0; i<size; i++)
758 fscanf(fd,"%lf %lf %d",&x,&y,&z);
763 SetKeyContours(&tempX,&tempY,&tempZ);
769 //---------------------------------------------------------------------------------------------------------
770 int PropContour::VectorDirection(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
773 double SumX = 0,SumY = 0;
774 double ax,ay,bx,by,axb;
775 int size = InX->size();
776 for(i=0; i<size; i++)
778 SumX = SumX + (*InX)[i];
779 SumY = SumY + (*InY)[i];
781 SumX = SumX/size; //Mass Center: X coord
782 SumY = SumY/size; //Mass Center: Y coord
786 for(i=0; i<size; i++)
790 bx = (*InX)[(i+1)%size]-SumX;
791 by = (*InY)[(i+1)%size]-SumY;
792 axb = (ax*by) - (bx*ay);
802 if(positive >= negative)
813 //----------------------------------------------------------------------------------------------------
814 void PropContour::VectorOrder(int dir, int posini, std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
816 int size = InX->size();
818 std::vector<double> tempX;
819 std::vector<double> tempY;
820 std::vector<double> tempZ;
826 for(i=0; i<size; i++)
830 tempX.push_back((*InX)[posini]);
831 tempY.push_back((*InY)[posini]);
832 tempZ.push_back((*InZ)[posini]);
841 tempX.push_back((*InX)[posini]);
842 tempY.push_back((*InY)[posini]);
843 tempZ.push_back((*InZ)[posini]);
854 for(i=0; i<size; i++)
856 InX->push_back(tempX[i]);
857 InY->push_back(tempY[i]);
858 InZ->push_back(tempZ[i]);
862 //----------------------------------------------------------------------------------------------------
863 void PropContour::PreparePointsForSpline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ,
864 std::vector<int>*Sizes)
866 int sizeS = Sizes->size();
867 //int sizeV = InX->size(); // JPRx
868 int i,j,mem,posinic,dir,cont;
871 std::vector<double> tempX;
872 std::vector<double> tempY;
873 std::vector<double> tempZ;
874 std::vector<double> lstX;
875 std::vector<double> lstY;
876 std::vector<double> lstZ;
884 for(i=0; i<sizeS; i++)
890 for(j=0; j<(*Sizes)[i]; j++)
892 tempX.push_back((*InX)[j+mem]);
893 tempY.push_back((*InY)[j+mem]);
894 tempZ.push_back((*InZ)[j+mem]);
895 if( (*InX)[j] < leX )
901 mem = mem + (*Sizes)[i];
902 dir = VectorDirection(&tempX,&tempY,&tempZ);
903 VectorOrder(dir,posinic,&tempX,&tempY,&tempZ);
905 for(j=0; j<(*Sizes)[i]; j++)
907 lstX.push_back(tempX[j]);
908 lstY.push_back(tempY[j]);
909 lstZ.push_back((*InZ)[j+cont]);
911 cont = cont + (*Sizes)[i];
914 //Fill the Finally lst in X,Y,Z ---------------
915 int sizetemp = lstX.size();
916 //printf("\nJSTG-PropContour::PreparePointsForSpline");
920 for(i=0; i<sizetemp; i++)
922 //printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]);
923 InX->push_back(lstX[i]);
924 InY->push_back(lstY[i]);
925 InZ->push_back(lstZ[i]);
930 //----------------------------------------------------------------------------------------------------
931 vtkImageData* PropContour::method_Spline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ, std::vector<int>*Sizes)
933 //long interval = wxGetElapsedTime(TRUE); // JPRx
935 int i,j,k,sizeX,sizeS,sizeInS;
938 //EED double spc[3]={1,1,1};
939 std::vector<double> lstX;
940 std::vector<double> lstY;
941 std::vector<double> lstZ;
942 std::vector<double> tempX;
943 std::vector<double> tempY;
944 std::vector<double> tempZ;
945 std::vector<double> tempS;
946 _mContourModel = new manualContourModel();
947 _mContourModel->SetNumberOfPointsSpline(_interpnumber);
951 // _dimImage[0] = 200; // X axis
952 // _dimImage[1] = 200; // Y axis
953 // _dimImage[2] = 200; // Z axis
955 // unsigned short *pValue;
956 // imagedataValue = vtkImageData::New();
957 // imagedataValue->SetScalarTypeToUnsignedShort();
958 // imagedataValue->SetSpacing(spc);
959 // imagedataValue->SetDimensions(_dimImage);
960 // imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
961 // imagedataValue->AllocateScalars();
962 // imagedataValue->Update();
969 for(i=0; i<sizeX; i++)
971 lstX.push_back((*InX)[i]);
972 lstY.push_back((*InY)[i]);
973 lstZ.push_back((*InZ)[i]);
976 PreparePointsForSpline(&lstX,&lstY,&lstZ,Sizes);
977 /*int sizetemp = lstX.size();
978 printf("\nJSTG-PropContour::method_Spline");
979 for(i=0; i<sizetemp; i++)
981 printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]);
984 sizeS = Sizes->size();
986 for(i=0; i<sizeS; i++)
988 _mContourModel->DeleteAllPoints();
989 _mContourModel->SetCloseContour(true);
990 sizeInS = (*Sizes)[i];
991 for(j=0; j<sizeInS; j++)
993 _mContourModel->AddPoint(lstX[cont],lstY[cont],lstZ[cont]);
996 _mContourModel->UpdateSpline();
997 numspline = _mContourModel->GetNumberOfPointsSpline();
998 for(k=0; k<numspline; k++)
1000 _mContourModel->GetSpline_i_Point(k,&x,&y,&z);
1001 //EED pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z);
1009 int tam = numspline;
1010 std::vector<Vector>::iterator it;
1014 for(i=0; i<numspline; i++)
1016 Vector *vec = new Vector();
1017 _planevector.push_back(*vec);
1020 for(j=0; j<tam; j++)
1022 _mContourModel->DeleteAllPoints();
1023 _mContourModel->SetCloseContour(false);
1025 for(i=0; i<sizeS; i++)
1027 //double hh=tempZ[cont+j]; // JPRx
1028 _mContourModel->AddPoint(tempX[cont+j],tempY[cont+j],tempZ[cont+j]);
1031 _mContourModel->UpdateSpline();
1032 numspline = _mContourModel->GetNumberOfPointsSpline();
1033 for(k=0; k<numspline; k++)
1036 _mContourModel->GetSpline_i_Point(k,&x,&y,&z);
1037 _planevector[k].set_x(x);
1038 _planevector[k].set_y(y);
1039 _planevector[k].set_z(z);
1040 _planevector[k].set_plane(k);
1041 //EED pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z);
1042 //EED *pValue = 128;
1047 int tempsize = _planevector.size();
1049 for(i=0; i<tempsize; i++)
1051 sizexx = _planevector[i].getsize_x();
1055 interval = wxGetElapsedTime(FALSE);
1056 int zplanes = abs(tempZ[0]-tempZ[tempZ.size()-1]);
1057 long interPlane = interval/zplanes;
1059 for(i=0; i<sizeS; i++)
1061 pointsSize = pointsSize + (*Sizes)[i];
1063 for(i=0; i<pointsSize; i++)
1065 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*InX)[i] , (int)(*InY)[i] ,(int)(*InZ)[i]);
1069 printf("\n\n JSTG - PropContour::method_Spline ------------------------");
1070 printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %d (ms)",interval);
1071 //printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
1072 printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %d (ms)",interPlane);
1073 printf("\n NUMBER OF PLANES................................. %d",zplanes);
1074 printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
1075 printf("\n ------------------------------------------------------------");
1077 return imagedataValue;
1080 //---------------------------------------------------------------------------------------------------
1081 void PropContour::ResetPlaneVector()
1084 // int ii,iiSize = _planevector.size();
1085 // for (ii=0 ; ii<iiSize ; ii++)
1087 // vec = &(_planevector[ii]);
1090 _planevector.clear();
1094 //---------------------------------------------------------------------------------------------------
1095 void PropContour::ResetKeyContours()
1097 _KeyContourSizes.clear();
1098 _KeyContourX.clear();
1099 _KeyContourY.clear();
1100 _KeyContourZ.clear();
1102 //---------------------------------------------------------------------------------------------------
1103 void PropContour::SetKeyContours(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
1106 int idKeyContour = 0;
1107 int idKeyContourSizes = 0;
1108 int tmpIdKeyContSizes = 0;
1109 bool okFind = false;
1111 int sizeKeyContour,Z=(int)(*InZ)[0];
1112 sizeKeyContour = _KeyContourZ.size();
1113 for (i=0; i<sizeKeyContour; i++)
1117 if ( (_KeyContourZ[i-1]<Z) && (_KeyContourZ[i]>=Z) )
1120 idKeyContourSizes = tmpIdKeyContSizes;
1124 if ( (i<sizeKeyContour) && (_KeyContourZ[i-1] != _KeyContourZ[i]) )
1126 tmpIdKeyContSizes++;
1129 if (_KeyContourZ[0]>Z)
1132 idKeyContourSizes = 0;
1141 idKeyContour = _KeyContourX.size();
1142 idKeyContourSizes = _KeyContourSizes.size();
1146 _KeyContourSizes.insert( _KeyContourSizes.begin() + idKeyContourSizes , InX->size() );
1147 for(i=0; i<(int)(InX->size()); i++)
1149 _KeyContourX.insert( _KeyContourX.begin() + idKeyContour, (*InX)[i] );
1150 _KeyContourY.insert( _KeyContourY.begin() + idKeyContour, (*InY)[i] );
1151 _KeyContourZ.insert( _KeyContourZ.begin() + idKeyContour, (*InZ)[i] );
1156 // _KeyContourSizes.push_back( InX->size() );
1157 // for(i=0; i<InX->size(); i++)
1159 // _KeyContourX.push_back( (*InX)[i] );
1160 // _KeyContourY.push_back( (*InY)[i] );
1161 // _KeyContourZ.push_back( (*InZ)[i] );
1165 //---------------------------------------------------------------------------------------------------
1166 vtkImageData* PropContour::CalculeSplinePropagation()
1168 if(_KeyContourSizes.size() <= 0)
1170 printf("\n There would be at last 1 contour");
1173 if(_KeyContourSizes.size() == 1)
1177 if(_KeyContourSizes.size() >= 2)
1179 return method_Spline(&_KeyContourX,&_KeyContourY,&_KeyContourZ,&_KeyContourSizes);
1183 //---------------------------------------------------------------------------------------------------
1184 void PropContour::GetKeyContours(std::vector<double>*KeyX, std::vector<double>*KeyY, std::vector<double>*KeyZ, std::vector<int>*KeyS)
1192 for(i=0; i<(int)(_KeyContourSizes.size()); i++)
1194 KeyS->push_back( _KeyContourSizes[i] );
1196 for(i=0; i<(int)(_KeyContourX.size()); i++)
1198 KeyX->push_back( _KeyContourX[i] );
1199 KeyY->push_back( _KeyContourY[i] );
1200 KeyZ->push_back( _KeyContourZ[i] );
1203 //---------------------------------------------------------------------------------------------------
1204 void PropContour::GetPropagatedContours( std::vector<Vector>*planevec )
1208 for(i=0; i<(int)(_planevector.size()); i++)
1210 Vector *temp = new Vector();
1211 temp->set_plane( _planevector[i].get_plane() );
1212 for(j=0; j<_planevector[i].getsize_x(); j++)
1214 temp->set_x( _planevector[i].get_x(j) );
1215 temp->set_y( _planevector[i].get_y(j) );
1216 temp->set_z( _planevector[i].get_z(j) );
1218 planevec->push_back(*temp);
1222 //---------------------------------------------------------------------------------------------------
1223 void PropContour::SetInterpNumber(int val)
1225 _interpnumber = val;
1228 int PropContour::FindIdWithZ(double z)
1231 int k,size=_planevector.size();
1232 double dist,minDist=99999999;
1233 for ( k=0 ; k<size ; k++)
1235 dist=fabs( z - _planevector[k].get_z(0) );
1245 //---------------------------------------------------------------------------------------------------
1246 void PropContour::GetIdContour(int id, std::vector<double>*vecX, std::vector<double>*vecY, std::vector<double>*vecZ)
1252 //int sizeplane = _planevector[id].getsize_x();
1254 for(i=0; i<_planevector[id].getsize_x(); i++)
1256 vecX->push_back( _planevector[id].get_x(i) );
1257 tempx = _planevector[id].get_x(i);
1258 vecY->push_back( _planevector[id].get_y(i) );
1259 vecZ->push_back( _planevector[id].get_z(i) );
1262 //---------------------------------------------------------------------------------------------------
1263 //---------------------------------------------------------------------------------------------------
1264 //---------------------------------------------------------------------------------------------------