vtkOtsuImageData.cxx

Go to the documentation of this file.

#include "vtkOtsuImageData.h"

         vtkOtsuImageData::vtkOtsuImageData(){
                 volume = NULL;
                 P = NULL;

        }

         vtkOtsuImageData::~vtkOtsuImageData(){
                
        }

        int vtkOtsuImageData::calculateOptimalThreshold(vtkImageData* voi){
                
                int extent[6];
                int static TAM = 10000;
                
                
                this->volume = voi; //VOI
                this->volume->GetExtent(extent);

                printf("otsu: extent = x (%d %d) y (%d %d) z(%d %d)\n",extent[0],extent[1],extent[2],extent[3],extent[4],extent[5]);

                                        
                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
                
                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++){
                                        
                                                                                
                                        valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
                                        if (valor > max){
                                                max = (int)valor;
                                        }
                                        
                                        index  = (int)valor;
                                        this->H[index]++;
                                        voxelCount++;
                                        
                                }
                        }
                }

                printf("otsu: voxelCount = %d\n",voxelCount);
                
                //3. Normalizacion de la funcion P
                for(index=0;index<=max;index++){
                        this->P[index] = (double)this->H[index] / (double)voxelCount;
                        //printf("H[%d] = %d\n",index, H[index]);
                }




                //FILE *logger = NULL;
        
                //if ((logger = freopen("c:\\Diego\\Andes\\Tesis\\ITK\\FM_3D\\bin\\Debug\\thotsu.txt","w", stdout)) == NULL)
                        //exit(-1);

                //printf("THRESHOLD; VARIANZA INTERCLASE\n");


                //4. Recorrer los diferentes niveles de intensidad

                for(T=0;T<=max;T++){

                        //printf("T= %d, ",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;
                        //printf("VI = %.2f\n",SI[T]);
                }
                
            int optimalThreshold        = 0;
                double optimalCriteria  = VTK_LARGE_FLOAT;

        

                for (int index2=0;index2<=max;index2++){
                        
                        if (SI[index2] < optimalCriteria){
                                optimalCriteria = SI[index2];
                                optimalThreshold = index2;
                                //printf("ENTRO %d %d\n",index2,optimalThreshold);
                        }
                        //printf("SI[%d] = %.2f, optimo = %d\n",index2,SI[index2],optimalThreshold);
                }
                //fclose(logger);
                //logger = NULL;

                printf(">>>%d<<<\n", optimalThreshold);
        
                return optimalThreshold;
        }