vtkOtsuSphereSource.cxx

Go to the documentation of this file.

#include "vtkOtsuSphereSource.h"

         vtkOtsuSphereSource::vtkOtsuSphereSource(){
                 volume = NULL;
                 region = NULL;
                 result = NULL;
                 P = NULL;

        }

         vtkOtsuSphereSource::~vtkOtsuSphereSource(){
                
        }

        double vtkOtsuSphereSource::calculateOptimalThreshold(vtkImageData* voi){
                
                int extent[6];
                double centro[3];
                int static TAM = 2000;

                        
                
                this->volume = voi; //VOI
                this->region = region; //JF
                this->volume->GetExtent(extent);

                centro[0] = (extent[0] + extent[1]) / 2;
                centro[1] = (extent[2] + extent[3]) / 2;
                centro[2] = (extent[4] + extent[5]) / 2;

                double radio = (extent[1] - extent[0]) /(double)2.0;

                result = vtkImageData::New();
        result->DeepCopy(volume);
                
                double valorEnRegion = 0.0;
                double valor = 0.0;
                int index = 0;
                int T = 0;
                int max=0;
                int voxelCount = 0;
        
                //1. Inicializacion
                this->H = new int[TAM];
                this->P = new double[TAM];
                this->SI = new double[TAM];
                for (index = 0; index <TAM; index++){
                        this->H[index] = 0;
                        this->P[index] = 0.0;   //probabilidad acumulada
                        this->SI[index] = 0.0;
                }
                
                //2. Calcular el Histograma
                double *ptr;
                int i, j, k;
                for(i=extent[0];i<=extent[1];i++){
                        for(j=extent[2];j<=extent[3];j++){
                                for(k=extent[4];k<=extent[5];k++){
                                        
                                        ptr = (double *)this->result->GetScalarPointer(i,j,k);

                                        // verificar si esta dentro de la esfera
                                        if(  ((i-centro[0])*(i-centro[0]) + 
                                                 (j-centro[1])*(j-centro[1]) +
                                                 (k-centro[2])*(k-centro[2])) <= (radio*radio)){

                                                
                                                        valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
                                                        if (valor > max){
                                                                max = (int)valor;
                                                        }
                                                        index  = (int)valor;
                                                        this->H[index]++;
                                                        voxelCount++;
                                                
                                                
                                        }
                                        else {
                                                *ptr=0.0;
                                        }
                                        
                                }
                        }
                }
                
                //3. Normalizacion de la funcion P
                for(index=0;index<=max;index++){
                        this->P[index] = (double)this->H[index] / (double)voxelCount;
                }


                //4. Recorrer los diferentes niveles de intensidad

                for(T=0;T<=max;T++){
                        
                        //4.1 Calcular las probabilidades acumuladas de las clases A y B
                        double Q1 = 0.0;
                        double Q2 = 0.0;
                        for(index=0;index<=max;index++){
                                if (index < T){
                                        Q1 = Q1 + this->P[index];
                                }
                                else{
                                        Q2 = Q2 + this->P[index];
                                }
                        }

                        //4.2 Calcular las medias de las clases A y B
                        double U1 = 0.0;
                        double U2 = 0.0;
                        for (index=0;index<=max;index++){
                                if (index < T){
                                        U1 = U1 + ((double)(index * this->P[index])/(double)Q1);
                                }
                                else{
                                        U2 = U2 + ((double)(index * this->P[index])/(double)Q2);
                                }
                        }

                        
                        //4.3 Calcular las varianzas individuales
                        double S1 = 0.0;
                        double S2 = 0.0;
                        for (index=0;index<=max;index++){
                                if (index < T){
                                        S1 = S1 + ((index - U1)*(index - U1)*(this->P[index]/(double)Q1));
                                }
                                else{
                                        S2 = S2 + ((index - U2)*(index - U2)*(this->P[index]/(double)Q2));
                                }
                        }

                        //4.4 Calcular SI: La varianza intraclases y seleccionar el valor minimo
                        this->SI[T] = Q1*S1 + Q2*S2;
                }
                
            int optimalThreshold        = 0;
                double optimalCriteria  = VTK_LARGE_FLOAT;

                for (index=0;index<=max;index++){
                        if (SI[index] < optimalCriteria){
                                optimalCriteria = SI[index];
                                optimalThreshold = index;
                        }
                }
                return optimalThreshold;
        }




        vtkImageData* vtkOtsuSphereSource::getImageForSegmentation(){
                return this->result;
        }