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;
        this->imagedata->GetScalarRange(scalarRange);

        // RaC 20-11-09 Changes in InitLstContoursLinesYPoints
        int ext[6];
        this->imagedata->GetWholeExtent(ext);
        _sizeImageY = ext[3]-ext[2]+1;

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

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


//------------------------------------------------------------------------
void ContourExtractData::GetMinMaxPoint(int *minPoint, 
                                                                                int *maxPoint, 
                                                                                manualBaseModel *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;
}



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

        /* RaC 20-11-09 (C1) Changes to use the method without the image.
        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++)
                {
                        // To process the statistics of the Points contour the procedure is different
                        // RaC 19-09-09
                        manualBaseModel *mbm = lstManConMod[i];
                        if(mbm->GetTypeModel()==7)
                        {
                                if(mbm->IsPoint(x,y)==true)
                                {
                                        numberLeft=1;
                                }
                        }//if
                        else
                        {
                                numberLeft =  AnalisisContourInsideV2(x,y, i );
                        }//else

                        if ( (numberLeft % 2) ==1){         numberInside++;  }
                }
                if ( numberInside == (size) ){ result=true; }
        } // AND  Intersection

        numberLeft=0;

        if (typeOperation==1)  // OR  All
        {
                for (i=0;i<size;i++)
                {
                        // To process the statistics of the Points contour the procedure is different
                        // RaC 19-09-09
                        manualBaseModel *mbm = lstManConMod[i];
                        if(mbm->GetTypeModel()==7)
                        {
                                if(mbm->IsPoint(x,y)==true)
                                {
                                        numberLeft=1;
                                }
                        }//if
                        else
                        {
                                numberLeft =  AnalisisContourInsideV2(x,y, i );
                        }//else
                        if ( (numberLeft % 2) ==1){ result=true;  }
                }
        } // OR  All

        numberLeft=0;

        if (typeOperation==2)  // XOR  crown
        {
                for (i=0;i<size;i++)
                {
                        // To process the statistics of the Points contour the procedure is different
                        // RaC 19-09-09
                        manualBaseModel *mbm = lstManConMod[i];
                        if(mbm->GetTypeModel()==7)
                        {
                                if(mbm->IsPoint(x,y)==true)
                                {
                                        numberLeft=1;
                                }
                        }//if
                        else
                        {
                                numberLeft =  numberLeft + AnalisisContourInsideV2(x,y, i );
                        }//else

                }
                if ( numberLeft % 2 ==1){       result = true;  } 
        }// XOR  crown



        //RaC }

        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++)
                        {

                                //RaC 20-11-09 Changes specified in isInside (C1)
                                int ext[6];
                                imagedata->GetExtent(ext);

                                if ((i>=0) && (i<=ext[1]) && (j>=0) && (j<=ext[3]))
                                {
                  if (isInside(i,j,_typeOperation)==true)
                                        {
                                          value = GetDataValue(i,j,zImage);
                                          if ( (value>=scalarRange[0]) && (value<=scalarRange[1]) )
                                          {
                                                PutVtkImageDataResultValue(i,j,zImage,  value );
                                          } // scalarRange
                                        } // if isInside

                                } // if ext
                        } // for i
                } // for j


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

}

//------------------------------------------------------------------------
void ContourExtractData::GetValuesInsideCrown(  int *numberOfPixels, 
                                                                                                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;
        int acum=0;

        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++)
                {
                        //RaC 20-11-09 Changes specified in isInside (C1)
                        int ext[6];
                        imagedata->GetExtent(ext);

                        if ((i>=0) && (i<=ext[1]) && (j>=0) && (j<=ext[3]))
                        {
                    if (isInside(i,j,_typeOperation)==true)
                                {

                                        acum++;
                                        value = GetDataValue(i,j,zImage);
                                        if ( (value>=scalarRange[0]) && (value<=scalarRange[1]) )
                                        {
                                                pLstValue               -> push_back( value );
                                                pLstValuePosX   -> push_back( i );
                                                pLstValuePosY   -> push_back( j );
                                                pLstValuePosZ   -> push_back( -1 );
                                        } // scalarRange
                                } // if isInside
                        } // if ext
                } // for
        } // for

        *numberOfPixels = acum;
}

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

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, (*rsize),
                                                                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;
        manualBaseModel *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();

        /* RaC 20-11-09
        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<_sizeImageY ; 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,_sizeImageY);
        }

}

void ContourExtractData::SetScalarRange(double min, double max)
{
        scalarRange[0]=min;
        scalarRange[1]=max;
}

void ContourExtractData::SetSizeImageY(int pSizeImageY)
{
        _sizeImageY=pSizeImageY;
}