2 #include "ContourExtractData.h"
5 //----------------------------------------------------------------------
6 ContourExtractData::ContourExtractData( bool okImagesResults)
8 this->imagedata = NULL;
9 imagedataValueResult = NULL;
10 imagedataMaskResult = NULL;
11 this->okImagesResults = okImagesResults;
15 // ------------------------------------------------------------------------
17 ContourExtractData::~ContourExtractData()
22 //----------------------------------------------------------------------
23 void ContourExtractData::SetImage( vtkImageData* imagedata)
25 this->imagedata = imagedata;
27 // init vtk image result : valuesImage maskImage
28 if (this->okImagesResults==true){ InitVtkImagesResult(); }
30 //----------------------------------------------------------------------
31 void ContourExtractData::SetZtoBeAnalys( int z )
36 //------------------------------------------------------------------------
37 void ContourExtractData::SetLstManualContourModel( std::vector<manualContourModel*> lstManConMod)
39 this->lstManConMod = lstManConMod;
43 //------------------------------------------------------------------------
44 void ContourExtractData::GetMinMaxPoint(int *minPoint,
46 manualContourModel *manualcontourmodel
50 //int np = manualcontourmodel->GetSizeLstPoints( ); // number of control points // JPRx
52 // JSTG 26-02-08 ---------------------------------------------------------------------------------------
53 //int nps = manualviewbaseecontour->GetNumberOfPointsSpline(); // number of points in the spline
54 int nps = manualcontourmodel->GetNumberOfPointsSpline(); // number of points in the spline
55 //------------------------------------------------------------------------------------------------------
57 // JSTG 26-02-08 ---------------------------------------------------------------------------------------
60 //double delta=( double ) ( np ) / ( double ) ( nps );
61 manualcontourmodel->UpdateSpline();
64 //t= delta * (double)i;
65 //manualcontourmodel->GetSplinePoint(t,x,y,z);
66 manualcontourmodel->GetSpline_i_Point(i,&x,&y,&z);
67 if (x<minPoint[0]){ minPoint[0]=(int)x; }
68 if (y<minPoint[1]){ minPoint[1]=(int)y; }
69 if (x>maxPoint[0]){ maxPoint[0]=(int)x; }
70 if (y>maxPoint[1]){ maxPoint[1]=(int)y; }
76 //------------------------------------------------------------------------------------------------------
79 //------------------------------------------------------------------------
80 void ContourExtractData::GetMinMaxPoint_Of_LstManConMod( int *minPoint,
85 int i,size = lstManConMod.size();
86 for(i=0 ; i<size ; i++)
88 GetMinMaxPoint(minPoint,maxPoint,lstManConMod[i]);
92 //------------------------------------------------------------------------
93 int ContourExtractData::AnalisisContourInside(int x,
95 manualContourModel *manualcontourmodel
101 //int np = manualcontourmodel->GetSizeLstPoints( ); // number of control points // JPRx
103 // JSTG 26-02-08 ---------------------------------------------------------------------------------------
104 //int nps = manualviewbaseecontour->GetNumberOfPointsSpline(); // number of points in the spline
105 int nps = manualcontourmodel->GetNumberOfPointsSpline(); // number of points in the spline
106 double x1,y1,z1,x2,y2,z2;
107 double borderX, borderY;
108 double xx1, yy1,xx2, yy2;
109 //double delta=( double ) ( np ) / ( double ) ( nps );
110 manualcontourmodel->UpdateSpline();
111 //------------------------------------------------------------------------------------------------------
116 // JSTG 26-02-08 ---------------------------------------------------------------------------------------
118 //manualcontourmodel->GetSplinePoint(0,x1,y1,z1);
119 manualcontourmodel->GetSpline_i_Point(0,&x1,&y1,&z1);
120 x1=x1+0.5; y1=y1+0.5;
121 for (i=1; i<=nps; i++)
127 //t= delta * (double)(i%nps);
128 //manualcontourmodel->GetSplinePoint(t,x2,y2,z2);
129 manualcontourmodel->GetSpline_i_Point(i,&x2,&y2,&z2);
130 x2=x2+0.5; y2=y2+0.5;
131 //------------------------------------------------------------------------------------------------------
133 //by triangle similarity
134 if ( ((y1<y2)&&(y>=y1)&&(y<y2)) || ((y1>y2)&&(y>=y2)&&(y<y1)) )
136 if (y1<y2) { xx1=x1; yy1=y1; xx2=x2; yy2=y2;} else { xx1=x2; yy1=y2; xx2=x1; yy2=y1; }
138 d = ( fabs(xx2-xx1)*(y-yy1) ) / (yy2-yy1) ;
139 if ( (xx1<xx2)&&(x<(xx1+d)) )
145 if ( (xx1>xx2)&&(x<(xx1-d)) ) {
150 } // if point inside contour
152 //Border verification
153 if ( (x==(int)borderX) && (y==(int)borderY) ) { inBorder=true; }// if point in border
155 // Verification : border in horizontal line
156 if ( ((int)y1==(int)y2) && ((int)y1==y) && (x1<x2) && (x>=x1) && (x<=x2)) { inBorder=true; }
157 if ( ((int)y1==(int)y2) && ((int)y1==y) && (x2<x1) && (x>=x2) && (x<=x1)) { inBorder=true; }
159 if (inBorder==true){ i=nps; }
166 if (inBorder==true) { result=1; }
172 //------------------------------------------------------------------------
173 int ContourExtractData::AnalisisContourInsideV2(int x, int y, int iContour )
179 int nps=_lstlstlstVecX1[iContour][y].size();
182 double borderX, borderY;
183 double xx1, yy1,xx2, yy2;
187 for (i=0; i<nps; i++)
189 x1=_lstlstlstVecX1[iContour][y][i];
190 y1=_lstlstlstVecY1[iContour][y][i];
191 x2=_lstlstlstVecX2[iContour][y][i];
192 y2=_lstlstlstVecY2[iContour][y][i];
199 xx1=x1; yy1=y1; xx2=x2; yy2=y2;
201 xx1=x2; yy1=y2; xx2=x1; yy2=y1;
204 double difxx2xx1=fabs(xx2-xx1);
205 if (difxx2xx1==0) difxx2xx1=0.0000000001;
207 // Finding border looking in vertical direction AND verifing if pixel is at right of the line
208 if ( (yy>=yy1)&&(yy<=yy2) )
211 //by triangle similarity
212 d = ( fabs(xx2-xx1)*(yy-yy1) ) / (yy2-yy1) ;
213 if ( (xx1<xx2)&&(x<(xx1+d)) ) { result++; }
214 if ( (xx1>xx2)&&(x<(xx1-d)) ) { result++; }
216 if ( (yy2-yy1)/difxx2xx1 >= 1.0)
227 } // if point inside y
230 // Finding border looking in vertical direction
231 if ( ((xx1<=xx2)&&(xx>=xx1)&&(xx<xx2)) || ((xx1>xx2)&&(xx>=xx2)&&(xx<xx1)) )
233 if ( (yy2-yy1)/difxx2xx1 <= 1.0)
235 //by triangle similarity
236 d = ( fabs(xx1-xx)*(yy2-yy1) ) / difxx2xx1;
242 } // if point inside x
245 //Border verification
246 if ( (x==(int)borderX) && (y==(int)borderY) ) { inBorder=true; }// if point in border
248 // Verification : border in horizontal line
249 if ( ((int)y1==(int)y2) && ((int)y1==y) && (x1<x2) && (x>=x1) && (x<=x2)) { inBorder=true; }
250 if ( ((int)y1==(int)y2) && ((int)y1==y) && (x2<x1) && (x>=x2) && (x<=x1)) { inBorder=true; }
252 if (inBorder==true){ i=nps; }
256 if (inBorder==true) { result=1; }
263 //------------------------------------------------------------------------
264 bool ContourExtractData::isInside(int x, int y, int typeOperation)
268 int i,size = this->lstManConMod.size();
269 int numberInside = 0;
272 imagedata->GetExtent(ext);
274 if ((x>=0) && (x<=ext[1]) && (y>=0) && (y<=ext[3]))
277 if (typeOperation==0) // AND Intersection
281 // numberLeft = AnalisisContourInside(x,y, lstManConMod[i] );
282 numberLeft = AnalisisContourInsideV2(x,y, i );
283 if ( (numberLeft % 2) ==1){ numberInside++; }
285 if ( numberInside == (size) ){ result=true; }
286 } // AND Intersection
289 if (typeOperation==1) // OR All
293 // numberLeft = AnalisisContourInside(x,y, lstManConMod[i] );
294 numberLeft = AnalisisContourInsideV2(x,y, i );
295 if ( (numberLeft % 2) ==1){ result=true; }
299 if (typeOperation==2) // XOR crown
303 // numberLeft = numberLeft + AnalisisContourInside(x,y, lstManConMod[i] );
304 numberLeft = numberLeft + AnalisisContourInsideV2(x,y, i );
306 if ( numberLeft % 2 ==1){ result = true; }
316 //------------------------------------------------------------------------
318 double ContourExtractData::GetDataValue(int x, int y, int z)
320 // wxVtk2DBaseView *wxvtk2dbaseview = (wxVtk2DBaseView*)wxvtkbaseview;
321 // int z = (int)wxvtk2dbaseview->GetVtkBaseData()->GetZ();
323 //EED OJO avec JS _zz = z;
327 p = imagedata->GetScalarPointer(x,y,z);
329 if (imagedata->GetScalarType()==VTK_CHAR)
332 result = (double)(*pp);
334 else if (imagedata->GetScalarType()==VTK_SIGNED_CHAR)
336 signed char *pp = (signed char*)p;
337 result = (double)(*pp);
339 else if (imagedata->GetScalarType()==VTK_UNSIGNED_CHAR)
341 unsigned char *pp = (unsigned char*)p;
342 result = (double)(*pp);
344 else if (imagedata->GetScalarType()==VTK_SHORT)
346 short *pp = (short*)p;
347 result = (double)(*pp);
349 else if (imagedata->GetScalarType()==VTK_UNSIGNED_SHORT)
351 unsigned short *pp = (unsigned short*)p;
352 result = (double)(*pp);
354 else if (imagedata->GetScalarType()==VTK_INT)
357 result = (double)(*pp);
359 else if (imagedata->GetScalarType()==VTK_UNSIGNED_INT)
361 unsigned int *pp = (unsigned int*)p;
362 result = (double)(*pp);
364 else if (imagedata->GetScalarType()==VTK_LONG)
367 result = (double)(*pp);
369 else if (imagedata->GetScalarType()==VTK_UNSIGNED_LONG)
371 unsigned long *pp = (unsigned long*)p;
372 result = (double)(*pp);
374 else if (imagedata->GetScalarType()==VTK_FLOAT)
376 float *pp = (float*)p;
377 result = (double)(*pp);
379 else if (imagedata->GetScalarType()==VTK_DOUBLE)
381 double *pp = (double*)p;
382 result = (double)(*pp);
388 //------------------------------------------------------------------------
390 void ContourExtractData::PutVtkImageDataResultValue( int x, int y, int z, double value )
392 unsigned short *pValue;
393 unsigned short *pMask;
394 pValue = (unsigned short *)imagedataValueResult->GetScalarPointer(x,y,z);
395 pMask = (unsigned short *)imagedataMaskResult->GetScalarPointer(x,y,z);
397 *pValue = (unsigned short)value;
400 //------------------------------------------------------------------------
401 void ContourExtractData::ResetImageResult(int z)
403 if (okImagesResults==true)
405 unsigned short *pValue;
406 unsigned short *pMask;
407 pValue = (unsigned short *)imagedataValueResult->GetScalarPointer(0,0,z);
408 pMask = (unsigned short *)imagedataMaskResult->GetScalarPointer(0,0,z);
411 imagedataValueResult->GetExtent(ext);
413 int size = (ext[1]-ext[0]+1) * (ext[3]-ext[2]+1);
414 memset(pValue,0,size*2);
415 memset(pMask,0,size*2);
420 //------------------------------------------------------------------------
421 void ContourExtractData::CalculateImageResult()
423 if (okImagesResults==true)
425 ResetImageResult(zImage);
432 minPoint[0] = 999999;
433 minPoint[1] = 999999;
434 maxPoint[0] = -999999;
435 maxPoint[1] = -999999;
437 GetMinMaxPoint_Of_LstManConMod(minPoint,maxPoint);
438 InitLstContoursLinesYPoints();
440 for (j=minPoint[1]; j<=maxPoint[1]; j++)
442 for (i=minPoint[0]; i<=maxPoint[0]; i++)
444 if (isInside(i,j,_typeOperation)==true)
446 value = GetDataValue(i,j,zImage);
447 PutVtkImageDataResultValue(i,j,zImage, value );
453 imagedataValueResult->Modified();
454 imagedataMaskResult->Modified();
455 imagedataValueResult->Update();
456 imagedataMaskResult->Update();
461 //------------------------------------------------------------------------
462 void ContourExtractData::GetValuesInsideCrown(std::vector<double> *pLstValue,
463 std::vector<double> *pLstValuePosX,
464 std::vector<double> *pLstValuePosY,
465 std::vector<double> *pLstValuePosZ)
468 pLstValuePosX->clear();
469 pLstValuePosY->clear();
470 pLstValuePosZ->clear();
472 // if (okImagesResults==true)
474 // ResetImageResult(zImage);
483 minPoint[0] = 999999;
484 minPoint[1] = 999999;
485 maxPoint[0] = -999999;
486 maxPoint[1] = -999999;
488 GetMinMaxPoint_Of_LstManConMod(minPoint,maxPoint);
489 InitLstContoursLinesYPoints();
491 for (j=minPoint[1]; j<=maxPoint[1]; j++)
493 for (i=minPoint[0]; i<=maxPoint[0]; i++)
495 if (isInside(i,j,_typeOperation)==true)
497 value = GetDataValue(i,j,zImage);
500 // if (okImagesResults==true){
501 // PutVtkImageDataResultValue(i,j,zImage, value );
504 pLstValue -> push_back( value );
505 pLstValuePosX -> push_back( i );
506 pLstValuePosY -> push_back( j );
507 pLstValuePosZ -> push_back( -1 );
514 // if (this->okImagesResults==true){
515 // imagedataValueResult->Modified();
516 // imagedataMaskResult->Modified();
522 //------------------------------------------------------------------------
524 vtkImageData *ContourExtractData::GetVtkImageValueResult()
526 return imagedataValueResult;
528 //------------------------------------------------------------------------
529 vtkImageData *ContourExtractData::GetVtkImageMaskResult()
531 return imagedataMaskResult;
533 // ------------------------------------------------------------------------
534 void ContourExtractData::InitVtkImagesResult()
541 imagedata->GetSpacing(spc);
542 imagedata->GetExtent(ext);
543 newDim[0]=ext[1]-ext[0]+1;
544 newDim[1]=ext[3]-ext[2]+1;
545 newDim[2]=ext[5]-ext[4]+1;
546 scalartype = imagedata->GetScalarType();
548 if (imagedataValueResult!=NULL)
550 imagedataValueResult->Delete();
552 imagedataValueResult = vtkImageData::New();
553 // imagedataValueResult->SetScalarType(scalartype);
554 imagedataValueResult->SetScalarTypeToUnsignedShort();
555 imagedataValueResult->SetSpacing(spc);
556 imagedataValueResult->SetDimensions( newDim );
557 imagedataValueResult->AllocateScalars();
559 if (imagedataMaskResult!=NULL)
561 imagedataMaskResult->Delete();
563 imagedataMaskResult = vtkImageData::New();
564 // imagedataMaskResult->SetScalarType(scalartype);
565 imagedataMaskResult->SetScalarTypeToUnsignedShort();
566 imagedataMaskResult->SetSpacing(spc);
567 imagedataMaskResult->SetDimensions( newDim );
568 imagedataMaskResult->AllocateScalars();
572 //------------------------------------------------------------------------
573 void ContourExtractData::InitVolumeStatistics()
577 vol_minValue = 9999999;
578 vol_maxValue =-9999999;
579 vol_acum_average = 0;
580 vol_acum_standardeviation = 0;
583 //------------------------------------------------------------------------
584 void ContourExtractData::SetVolumeStatistics(int rCountRange,
589 double acum_standardeviation)
591 vol_rCountRange = vol_rCountRange + rCountRange;
592 vol_rsize = vol_rsize + rsize;
594 if (minValue<vol_minValue){ vol_minValue = minValue; }
595 if (maxValue>vol_maxValue){ vol_maxValue = maxValue; }
597 vol_acum_average = vol_acum_average + acum_average;
598 vol_acum_standardeviation = vol_acum_standardeviation + acum_standardeviation;
601 //------------------------------------------------------------------------
602 void ContourExtractData::GetVolumeStatistics(int *vol_rCountRange,
604 double *vol_minValue,
605 double *vol_maxValue,
607 double *vol_standardeviation)
609 *vol_rCountRange = this->vol_rCountRange;
610 *vol_rsize = this->vol_rsize;
611 *vol_minValue = this->vol_minValue;
612 *vol_maxValue = this->vol_maxValue;
613 *vol_average = this->vol_acum_average / this->vol_rsize;
614 *vol_standardeviation = sqrt(this->vol_acum_standardeviation / this->vol_rsize);
618 //------------------------------------------------------------------------
619 void ContourExtractData::Statistics( std::vector<double> *inputLstValue,
627 double *rstandardeviation
633 double standardeviation = 0;
634 double acum_average = 0;
635 double acum_standardeviation = 0;
640 if (inputLstValue!=NULL)
642 size=inputLstValue->size();
644 max=(*inputLstValue)[0];
645 min=(*inputLstValue)[0];
646 // Average , countRange
648 for ( i=0; i<size; i++ )
650 ng=(*inputLstValue)[i];
651 acum_average = acum_average + ng;
652 if (max<ng) max=ng; // Max
653 if (min>ng) min=ng; // Min
654 if ((ng>=grayRangeMin) && (ng<=grayRangeMax)) countRange++; // countRange
656 average = acum_average / size;
659 acum_standardeviation=0;
661 for ( i=0; i<size; i++ )
663 tmp = (*inputLstValue)[i] - average;
664 acum_standardeviation = acum_standardeviation + tmp*tmp;
665 } // for standar deviation
666 standardeviation = sqrt(acum_standardeviation/size);
667 SetVolumeStatistics(countRange, size,
669 acum_average,acum_standardeviation);
675 *rCountRange = countRange;
679 *rstandardeviation = standardeviation;
682 //------------------------------------------------------------------------
683 void ContourExtractData::SetTypeOperation(int type)
688 //------------------------------------------------------------------------
689 void ContourExtractData::Fill_lstlstlstVecXY(int iContour, int sizeY)
692 double x1,y1,z1,x2,y2,z2;
693 manualContourModel *manualcontourmodel= lstManConMod[iContour];
694 int nps = manualcontourmodel->GetNumberOfPointsSpline(); // number of points in the spline
695 manualcontourmodel->UpdateSpline();
696 //------------------------------------------------------------------------------------------------------
698 for (y=0;y<sizeY;y++)
700 manualcontourmodel->GetSpline_i_Point(0,&x1,&y1,&z1);
701 x1=x1+0.5; y1=y1+0.5;
702 for (i=1; i<=nps; i++)
704 manualcontourmodel->GetSpline_i_Point(i,&x2,&y2,&z2);
705 x2=x2+0.5; y2=y2+0.5;
706 if ( ((y1<y2)&&(y>=y1)&&(y<=y2)) || ((y1>y2)&&(y>=y2)&&(y<=y1)) || ((int)y1==y) || ((int)y2==y) )
708 _lstlstlstVecX1[iContour][y].push_back(x1);
709 _lstlstlstVecY1[iContour][y].push_back(y1);
710 _lstlstlstVecX2[iContour][y].push_back(x2);
711 _lstlstlstVecY2[iContour][y].push_back(y2);
714 } // for i Points in spline
719 void ContourExtractData::InitLstContoursLinesYPoints()
721 // init InInside Optimisation
724 _lstlstlstVecX1.clear();
725 _lstlstlstVecY1.clear();
726 _lstlstlstVecX2.clear();
727 _lstlstlstVecY2.clear();
730 this->imagedata->GetWholeExtent(ext);
731 int sizeY = ext[3]-ext[2]+1;
732 std::vector<double> vecDouble;
733 std::vector< std::vector<double> > vecVecDouble;
734 for ( i=0 ; i<sizeY ; i++ )
736 vecVecDouble.push_back( vecDouble );
739 //Fill structure with points
740 int sizeContours = lstManConMod.size();
741 for( i=0 ; i<sizeContours ; i++ )
743 _lstlstlstVecX1.push_back( vecVecDouble );
744 _lstlstlstVecY1.push_back( vecVecDouble );
745 _lstlstlstVecX2.push_back( vecVecDouble );
746 _lstlstlstVecY2.push_back( vecVecDouble );
747 Fill_lstlstlstVecXY(i,sizeY);