ContourExtractData.cxx

Go to the documentation of this file.

#include "ContourExtractData.h"


//----------------------------------------------------------------------
ContourExtractData::ContourExtractData( bool okImagesResults)
  {
    this->imagedata                     = NULL;
        imagedataValueResult    = NULL;
        imagedataMaskResult             = NULL;
        this->okImagesResults   = okImagesResults;
        _typeOperation                  = 0;
  }

  // ------------------------------------------------------------------------

  ContourExtractData::~ContourExtractData()
  {
  }


//----------------------------------------------------------------------
void ContourExtractData::SetImage( vtkImageData* imagedata)
  {
        this->imagedata                 = imagedata;
                
        // init vtk image result : valuesImage maskImage  
        if (this->okImagesResults==true){  InitVtkImagesResult(); }
  }
//----------------------------------------------------------------------
void ContourExtractData::SetZtoBeAnalys( int z )
  {
        this->zImage                    = z;
  }

//------------------------------------------------------------------------
  void ContourExtractData::SetLstManualContourModel( std::vector<manualContourModel*> lstManConMod)
  {
          this->lstManConMod = lstManConMod;
  }


//------------------------------------------------------------------------
void ContourExtractData::GetMinMaxPoint(int *minPoint, 
                                                                                  int *maxPoint, 
                                                                                  manualContourModel *manualcontourmodel
                                                                                  )
{
        int i;
        //int   np              = manualcontourmodel->GetSizeLstPoints( );  // number of control points // JPRx

// JSTG 26-02-08 ---------------------------------------------------------------------------------------
        //int nps = manualviewbaseecontour->GetNumberOfPointsSpline(); // number of points in the spline
        int nps = manualcontourmodel->GetNumberOfPointsSpline(); // number of points in the spline
//------------------------------------------------------------------------------------------------------
        
// JSTG 26-02-08 ---------------------------------------------------------------------------------------
        //double x,y,z,t;
        double x,y,z;
        //double delta=( double ) ( np  ) / ( double ) ( nps  );
        manualcontourmodel->UpdateSpline();
        for (i=0; i<nps; i++)
        {
                //t= delta * (double)i;
                //manualcontourmodel->GetSplinePoint(t,x,y,z);
                manualcontourmodel->GetSpline_i_Point(i,&x,&y,&z);
                if (x<minPoint[0]){ minPoint[0]=(int)x; }
                if (y<minPoint[1]){ minPoint[1]=(int)y; }
                if (x>maxPoint[0]){ maxPoint[0]=(int)x; }
                if (y>maxPoint[1]){ maxPoint[1]=(int)y; }
        }
        minPoint[0]--;
        minPoint[1]--;
        maxPoint[0]++;
        maxPoint[1]++;
//------------------------------------------------------------------------------------------------------
}

//------------------------------------------------------------------------
void ContourExtractData::GetMinMaxPoint_Of_LstManConMod(        int *minPoint, 
                                                                                        int *maxPoint
                                                                                )

{
        int i,size = lstManConMod.size();
        for(i=0 ; i<size ; i++)
        {
                GetMinMaxPoint(minPoint,maxPoint,lstManConMod[i]);
        }
}


//------------------------------------------------------------------------
int ContourExtractData::AnalisisContourInsideV2(int x, int y, int iContour )
{
        bool inBorder=false;
        int result      = 0;
        int i;

        int nps=_lstlstlstVecX1[iContour][y].size();
        
        double x1,y1,x2,y2;
        double borderX, borderY;
        double xx1, yy1,xx2, yy2;
        double xx=x, yy=y;
        double d;

        for (i=0; i<nps; i++)
        {
                x1=_lstlstlstVecX1[iContour][y][i];
                y1=_lstlstlstVecY1[iContour][y][i];
                x2=_lstlstlstVecX2[iContour][y][i];
                y2=_lstlstlstVecY2[iContour][y][i];

                borderX=x1;
                borderY=y1;
                
                if (y1<y2) 
                { 
                        xx1=x1; yy1=y1; xx2=x2; yy2=y2;
                } else {
                        xx1=x2; yy1=y2; xx2=x1; yy2=y1; 
                } 
                
                double difxx2xx1=fabs(xx2-xx1);
                if (difxx2xx1==0) difxx2xx1=0.0000000001;
                
                // Finding border looking in vertical direction  AND  verifing if pixel is at right of the line 
                if  ( (yy>=yy1)&&(yy<=yy2) ) 
                {
                        //by triangle similarity
                        d = ( fabs(xx2-xx1)*(yy-yy1) ) / (yy2-yy1) ;
                        if ( (xx1<=xx2)&&(x<(xx1+d)) )  {       result++;       }
                        if ( (xx1>xx2)&&(x<(xx1-d)) )   {       result++;       } 
                        
                        if ( (yy2-yy1)/difxx2xx1 >= 1.0)
                        {
                                if (xx1<=xx2)   
                                { 
                                        borderX = xx1+d;
                                        borderY = y;
                                } else { 
                                        borderX = xx1-d;
                                        borderY = y;
                                } 
                        } 
                } // if point inside y 
                
                
                // Finding border looking in vertical direction
                if ( ((xx1<=xx2)&&(xx>=xx1)&&(xx<xx2)) || ((xx1>xx2)&&(xx>=xx2)&&(xx<xx1)) )
                {
                        if ( (yy2-yy1)/difxx2xx1 <= 1.0)
                        {
                                //by triangle similarity
                                d = ( fabs(xx1-xx)*(yy2-yy1) ) / difxx2xx1;
                                if (yy1+d<=yy2) {
                                        borderX=x;
                                        borderY=yy1+d;
                                }
                        }                       
                } // if point inside x                  
                
                
                //Border verification
                if (   (x==(int)borderX) && (y==(int)borderY)   )  { inBorder=true; }// if point in border
                
                // Verification : border in horizontal line 
                if ( ((int)y1==(int)y2)  &&  ((int)y1==y) && (x1<x2) && (x>=x1) && (x<=x2))             { inBorder=true;        }
                if ( ((int)y1==(int)y2)  &&  ((int)y1==y) && (x2<x1) && (x>=x2) && (x<=x1))             { inBorder=true;        }
                
                if (inBorder==true){ i=nps; }           
        } // for i
        
        if (inBorder==true)     { result=1;     }
        
        return result;
}



//------------------------------------------------------------------------
bool ContourExtractData::isInside(int x, int y, int typeOperation)
{
        bool result                             = false;
        int numberLeft                  = 0;
        int i,size                              = this->lstManConMod.size();
        int numberInside        = 0;

        int ext[6];
        imagedata->GetExtent(ext);

        if ((x>=0) && (x<=ext[1]) && (y>=0) && (y<=ext[3]))
        {

                if (typeOperation==0)  // AND  Intersection
                {
                        for (i=0;i<size;i++)
                        {
                                numberLeft =  AnalisisContourInsideV2(x,y, i );
                                if ( (numberLeft % 2) ==1){         numberInside++;  }
                        }
                        if ( numberInside == (size) ){ result=true; }
                } // AND  Intersection


                if (typeOperation==1)  // OR  All
                {
                        for (i=0;i<size;i++)
                        {
                                numberLeft =  AnalisisContourInsideV2(x,y, i );
                                if ( (numberLeft % 2) ==1){ result=true;  }
                        }
                } // OR  All

                if (typeOperation==2)  // XOR  crown
                {
                        for (i=0;i<size;i++)
                        {
                                numberLeft = numberLeft + AnalisisContourInsideV2(x,y, i );
                        }
                        if ( numberLeft % 2 ==1){       result = true;  } 
                }// XOR  crown



        }

        return result;
}

//------------------------------------------------------------------------

double ContourExtractData::GetDataValue(int x, int y, int z)
{
//      wxVtk2DBaseView *wxvtk2dbaseview = (wxVtk2DBaseView*)wxvtkbaseview;
//      int z = (int)wxvtk2dbaseview->GetVtkBaseData()->GetZ();
        //JSTG 13-03-08-----
//EED OJO avec JS       _zz = z;
        //------------------
        double result;
        void *p;
        p = imagedata->GetScalarPointer(x,y,z);

        if (imagedata->GetScalarType()==VTK_CHAR)
        {
                char *pp = (char*)p;
                result = (double)(*pp);
        }
        else if (imagedata->GetScalarType()==VTK_SIGNED_CHAR)
        {
                signed char *pp = (signed char*)p;
                result = (double)(*pp);
        }
        else if (imagedata->GetScalarType()==VTK_UNSIGNED_CHAR)
        {
                unsigned char *pp = (unsigned char*)p;
                result = (double)(*pp);
        }
        else if (imagedata->GetScalarType()==VTK_SHORT)
        {
                short *pp = (short*)p;
                result = (double)(*pp);
        }
        else if (imagedata->GetScalarType()==VTK_UNSIGNED_SHORT)
        {
                unsigned short *pp = (unsigned short*)p;
                result = (double)(*pp);
        }
        else if (imagedata->GetScalarType()==VTK_INT)
        {
                int *pp = (int*)p;
                result = (double)(*pp);
        }
        else if (imagedata->GetScalarType()==VTK_UNSIGNED_INT)
        {
                unsigned int *pp = (unsigned int*)p;
                result = (double)(*pp);
        }
        else if (imagedata->GetScalarType()==VTK_LONG)
        {
                long *pp = (long*)p;
                result = (double)(*pp);
        }
        else if (imagedata->GetScalarType()==VTK_UNSIGNED_LONG)
        {
                unsigned long *pp = (unsigned long*)p;
                result = (double)(*pp);
        }
        else if (imagedata->GetScalarType()==VTK_FLOAT)
        {       
                float *pp = (float*)p;
                result = (double)(*pp);
        }
        else if (imagedata->GetScalarType()==VTK_DOUBLE)
        {
                double *pp = (double*)p;
                result = (double)(*pp);
        }

        return result;
}

//------------------------------------------------------------------------

void ContourExtractData::PutVtkImageDataResultValue( int x, int y, int z, double value )
{
        unsigned short *pValue;
        unsigned short *pMask;
        pValue  = (unsigned short *)imagedataValueResult->GetScalarPointer(x,y,z);
        pMask   = (unsigned short *)imagedataMaskResult->GetScalarPointer(x,y,z);
        *pMask  = 255;
        *pValue = (unsigned short)value;
}

//------------------------------------------------------------------------
void ContourExtractData::ResetImageResult(int z)
{
        if (okImagesResults==true)
        {
                unsigned short *pValue;
                unsigned short *pMask;
                pValue  = (unsigned short *)imagedataValueResult->GetScalarPointer(0,0,z);
                pMask   = (unsigned short *)imagedataMaskResult->GetScalarPointer(0,0,z);
                
                int ext[6];
                imagedataValueResult->GetExtent(ext);

                int size = (ext[1]-ext[0]+1) * (ext[3]-ext[2]+1); 
                memset(pValue,0,size*2);
                memset(pMask,0,size*2);
        } // if
}


//------------------------------------------------------------------------
void ContourExtractData::CalculateImageResult()
{
        if (okImagesResults==true)
        {
                ResetImageResult(zImage);

                int minPoint[2];
                int maxPoint[2];
                int i,j;
                double value;

                minPoint[0] = 999999;
                minPoint[1] = 999999;
                maxPoint[0] = -999999;
                maxPoint[1] = -999999;

                GetMinMaxPoint_Of_LstManConMod(minPoint,maxPoint);
                InitLstContoursLinesYPoints();
                
                for (j=minPoint[1]; j<=maxPoint[1]; j++)
                {
                        for (i=minPoint[0]; i<=maxPoint[0]; i++)
                        {
                                if (isInside(i,j,_typeOperation)==true)
                                {
                                        value = GetDataValue(i,j,zImage);
                                        PutVtkImageDataResultValue(i,j,zImage,  value );
                                } // if
                        } // for i
                } // for j


                imagedataValueResult->Modified();
                imagedataMaskResult->Modified();
                imagedataValueResult->Update();
                imagedataMaskResult->Update();
        } // if

}

//------------------------------------------------------------------------
void ContourExtractData::GetValuesInsideCrown(std::vector<double> *pLstValue,
                                                                                                std::vector<double> *pLstValuePosX,
                                                                                                std::vector<double> *pLstValuePosY,
                                                                                                std::vector<double> *pLstValuePosZ)
{
        pLstValue->clear();
        pLstValuePosX->clear();
        pLstValuePosY->clear();
        pLstValuePosZ->clear();

//      if (okImagesResults==true)
//      {
//              ResetImageResult(zImage);
//      }

        int minPoint[2];
        int maxPoint[2];
        int i,j;
        double value;


        minPoint[0] = 999999;
        minPoint[1] = 999999;
        maxPoint[0] = -999999;
        maxPoint[1] = -999999;

        GetMinMaxPoint_Of_LstManConMod(minPoint,maxPoint);
    InitLstContoursLinesYPoints();

        for (j=minPoint[1]; j<=maxPoint[1]; j++)
        {
                for (i=minPoint[0]; i<=maxPoint[0]; i++)
                {
                        if (isInside(i,j,_typeOperation)==true)
                        {
                                value = GetDataValue(i,j,zImage);

// Borrame
//                              if (okImagesResults==true){
//                                      PutVtkImageDataResultValue(i,j,zImage,  value );
//                              }

                                pLstValue               -> push_back( value );
                                pLstValuePosX   -> push_back( i );
                                pLstValuePosY   -> push_back( j );
                                pLstValuePosZ   -> push_back( -1 );
                        } // if
                } // for
        } // for


// Borrame
//      if (this->okImagesResults==true){
//              imagedataValueResult->Modified();
//              imagedataMaskResult->Modified();
//      }


}

//------------------------------------------------------------------------

vtkImageData *ContourExtractData::GetVtkImageValueResult()
{
        return imagedataValueResult;
}
//------------------------------------------------------------------------
vtkImageData *ContourExtractData::GetVtkImageMaskResult()
{
        return imagedataMaskResult;
}
// ------------------------------------------------------------------------
void ContourExtractData::InitVtkImagesResult()
{
          int ext[6];
          int newDim[3];
          double spc[3];
          int scalartype;

          imagedata->GetSpacing(spc);
          imagedata->GetExtent(ext);
          newDim[0]=ext[1]-ext[0]+1;
          newDim[1]=ext[3]-ext[2]+1;
          newDim[2]=ext[5]-ext[4]+1;
          scalartype = imagedata->GetScalarType();

          if (imagedataValueResult!=NULL)
          {
                  imagedataValueResult->Delete();
          }
          imagedataValueResult = vtkImageData::New();
//        imagedataValueResult->SetScalarType(scalartype);
          imagedataValueResult->SetScalarTypeToUnsignedShort();
          imagedataValueResult->SetSpacing(spc);
          imagedataValueResult->SetDimensions( newDim );
          imagedataValueResult->AllocateScalars();

          if (imagedataMaskResult!=NULL)
          {
                  imagedataMaskResult->Delete();
          }
          imagedataMaskResult  = vtkImageData::New();
//        imagedataMaskResult->SetScalarType(scalartype);
          imagedataMaskResult->SetScalarTypeToUnsignedShort();
          imagedataMaskResult->SetSpacing(spc);
          imagedataMaskResult->SetDimensions( newDim );
          imagedataMaskResult->AllocateScalars();
}


//------------------------------------------------------------------------
void ContourExtractData::InitVolumeStatistics()
{
        vol_rCountRange                         = 0;
        vol_rsize                                       = 0;
        vol_minValue                            = 9999999;
        vol_maxValue                            =-9999999;
        vol_acum_average                        = 0;
        vol_acum_standardeviation       = 0;
}

//------------------------------------------------------------------------
void ContourExtractData::SetVolumeStatistics(int rCountRange, 
                                                                                        int rsize,
                                                                                        double minValue,
                                                                                        double maxValue,
                                                                                        double acum_average,
                                                                                        double acum_standardeviation)
{
        vol_rCountRange                         = vol_rCountRange + rCountRange; 
        vol_rsize                                       = vol_rsize + rsize; 
        
        if (minValue<vol_minValue){ vol_minValue = minValue;  }
        if (maxValue>vol_maxValue){ vol_maxValue = maxValue;  }
        
        vol_acum_average                        = vol_acum_average + acum_average; 
        vol_acum_standardeviation       = vol_acum_standardeviation + acum_standardeviation; 
}

//------------------------------------------------------------------------
void ContourExtractData::GetVolumeStatistics(int *vol_rCountRange, 
                                                                                         int *vol_rsize,
                                                                                         double *vol_minValue,
                                                                                         double *vol_maxValue,
                                                                                         double *vol_average,
                                                                                         double *vol_standardeviation)
{
        *vol_rCountRange                = this->vol_rCountRange; 
        *vol_rsize                              = this->vol_rsize; 
        *vol_minValue                   = this->vol_minValue; 
        *vol_maxValue                   = this->vol_maxValue; 
        *vol_average                    = this->vol_acum_average / this->vol_rsize;
        *vol_standardeviation   = sqrt(this->vol_acum_standardeviation / this->vol_rsize);
}


//------------------------------------------------------------------------
void ContourExtractData::Statistics( std::vector<double> *inputLstValue, 
                                                                        int     grayRangeMin,
                                                                        int     grayRangeMax,
                                                                        int             *rCountRange, 
                                                                        int             *rsize, 
                                                                        double  *rmin, 
                                                                        double  *rmax,
                                                                        double  *raverage,
                                                                        double  *rstandardeviation
                                                                                )
{
          double min                                            = 0;
          double max                                            = 0;
          double average                                        = 0;
          double standardeviation                       = 0;
          double acum_average                           = 0;
          double acum_standardeviation          = 0;
          int    size                                           = 0;
          int    countRange                                     = 0;
          double ng;

          if (inputLstValue!=NULL)
          {
                size=inputLstValue->size();
                if (size>0){
               max=(*inputLstValue)[0];
                   min=(*inputLstValue)[0];
     // Average , countRange
                        int i;
                        for ( i=0; i<size; i++ )
                        {
                                ng=(*inputLstValue)[i];
                                acum_average = acum_average + ng;
                                if (max<ng) max=ng;             // Max
                                if (min>ng) min=ng;     // Min
                                if ((ng>=grayRangeMin) && (ng<=grayRangeMax)) countRange++;  // countRange
                        } // for average
                        average = acum_average / size;

          // Standar Deviation
                        acum_standardeviation=0;
                        double tmp;
                        for ( i=0; i<size; i++ )
                        {
                tmp = (*inputLstValue)[i] - average;
                                acum_standardeviation = acum_standardeviation + tmp*tmp;
                        } // for standar deviation
                        standardeviation = sqrt(acum_standardeviation/size);
                        SetVolumeStatistics(countRange, size,
                                                                min,max,
                                                            acum_average,acum_standardeviation);
                } // if size
          } // if NULL

          // OUTPUT
                *rsize                          = size; 
                *rCountRange            = countRange;
                *rmin                           = min; 
                *rmax                           = max;
                *raverage                       = average;
                *rstandardeviation      = standardeviation;
}

//------------------------------------------------------------------------
void ContourExtractData::SetTypeOperation(int type)
{
        _typeOperation=type;
}

//------------------------------------------------------------------------
void ContourExtractData::Fill_lstlstlstVecXY(int iContour, int sizeY)
{
        int i,y;
        double x1,y1,z1,x2,y2,z2;
        manualContourModel *manualcontourmodel= lstManConMod[iContour]; 
        int nps = manualcontourmodel->GetNumberOfPointsSpline(); // number of points in the spline
        manualcontourmodel->UpdateSpline();
        //------------------------------------------------------------------------------------------------------
        
        for (y=0;y<sizeY;y++)
        {
                manualcontourmodel->GetSpline_i_Point(0,&x1,&y1,&z1);
                x1=x1+0.5; y1=y1+0.5;
                for (i=1; i<nps; i++)
                {
                        manualcontourmodel->GetSpline_i_Point(i,&x2,&y2,&z2);
                        x2=x2+0.5; y2=y2+0.5;
                        if ( ((y1<y2)&&(y>=y1)&&(y<=y2)) || ((y1>y2)&&(y>=y2)&&(y<=y1)) || ((int)y1==y) || ((int)y2==y)  )
                        {
                                _lstlstlstVecX1[iContour][y].push_back(x1);
                                _lstlstlstVecY1[iContour][y].push_back(y1);
                                _lstlstlstVecX2[iContour][y].push_back(x2);
                                _lstlstlstVecY2[iContour][y].push_back(y2);
                        } 
                        x1=x2; y1=y2; z1=z2;
                } // for i  Points in spline
        } // y

}

void ContourExtractData::InitLstContoursLinesYPoints()
{
        // init InInside Optimisation  
        int i;
        
        _lstlstlstVecX1.clear();
        _lstlstlstVecY1.clear();
        _lstlstlstVecX2.clear();
        _lstlstlstVecY2.clear();
        
        int ext[6];
        this->imagedata->GetWholeExtent(ext);
        int sizeY = ext[3]-ext[2]+1;
        std::vector<double> vecDouble;
        std::vector< std::vector<double> > vecVecDouble;
        for ( i=0 ; i<sizeY ; i++ )
        {
                vecVecDouble.push_back( vecDouble );
        }
        
        //Fill structure with points
        int sizeContours = lstManConMod.size();
        for( i=0 ; i<sizeContours ; i++ )
        {
                _lstlstlstVecX1.push_back( vecVecDouble );
                _lstlstlstVecY1.push_back( vecVecDouble );
                _lstlstlstVecX2.push_back( vecVecDouble );
                _lstlstlstVecY2.push_back( vecVecDouble );
                Fill_lstlstlstVecXY(i,sizeY);
        }
        
}

---