vtkOtsu.cxx

Go to the documentation of this file.

#include "vtkOtsu.h"

         vtkOtsu::vtkOtsu(){
                 volume = NULL;
                 region = NULL;
        
        }

         vtkOtsu::~vtkOtsu(){
                
        }




        double vtkOtsu::calculateOptimalThreshold(vtkImageData* volume, vtkImageData *region){

                this->volume = volume;
                this->region = region;
        
                double valor = 0.0;
                double media = 0.0;
                double valorEnRegion = 0.0;
            
        
                int extent[6];
                
                int index = 0;
                
                
                this->volume->GetExtent(extent);
                this->max=0.0;
                this->sigmaT = 0.0;
                
                //1. Numero de voxeles
                this->voxelCount = 0;
                


                //2. Inicializar el histograma y criterios
                this->histogram = new int[1000];
                this->criterias = new double[1000];
                for (index = 0; index <1000; index++){
                        this->histogram[index] = 0;
                        this->criterias[index] = 0.0;
                }
                
                
                //3. Encontrar el maximo, el promedio y 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++){
                                        valorEnRegion = this->region->GetScalarComponentAsDouble(i,j,k,0);
                                        if (valorEnRegion == 0.0){
                                                valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
                                                if (valor > this->max){
                                                        this->max = valor;
                                                }
                                                index  = (int)valor;
                                                this->histogram[index]++;
                                                media += valor;
                                                this->voxelCount++;
                                        }
                                }
                        }
                }
                
                media = (double)media/(double)this->voxelCount;

                //calcula sigma total
                for(index=0;index<=max;index++){
                        this->sigmaT += (index - media)*(index - media)* this->histogram[index] / (double)this->voxelCount;
                }


                //calcular criterio para cada posible threshold e ir seleccionando el maximo
            int optimalThreshold = 0;
                double optimalCriteria = 0.0;

                for (index=0;index<=max;index++){
                        criterias[index] = this->getCriteria(index);
                        if (criterias[index] > optimalCriteria){
                                optimalCriteria = criterias[index];
                                optimalThreshold = index;
                        }
                }
                return optimalThreshold;

                
        }


        double vtkOtsu::getCriteria(double threshold){
                
                int extent[6];
                int i,j,k;
                double valor = 0.0;
                double valorEnRegion = 0.0;

                this->volume->GetExtent(extent);
                
                //1. Encontrar la frecuencia y la media de las clases A y B
                double mediaA=0.0;
                double wA = 0.0;
                int countA = 0;
                
                double mediaB=0.0;
                double wB = 0.0;
                int countB = 0;

                for(i=extent[0];i<=extent[1];i++){
                        for(j=extent[2];j<=extent[3];j++){
                                for(k=extent[4];k<=extent[5];k++){
                                        valorEnRegion = this->region->GetScalarComponentAsDouble(i,j,k,0);
                                        if (valorEnRegion == 0.0){
                                                valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
                                                if (valor < threshold){
                                                        mediaA += valor;
                                                        countA++;

                                                }
                                                else{
                                                        mediaB += valor;
                                                        countB++;
                                                }
                                        }
                                }
                        }
                }

                mediaA = (double)mediaA / (double)countA;
                wA = (double)countA / (double)this->voxelCount;

                mediaB = (double)mediaB / (double)countB;
                wB = (double)countB / (double)this->voxelCount;


                //2. Encontrar la varianza de las clases A y B

                double sigmaA=0.0;
                double sigmaB=0.0;

                int index = 0;
                
                for(index=0;index<=threshold;index++){
                        sigmaA += (index - mediaA)*(index - mediaA)*this->histogram[index]/(double)countA;
                }

                
                for(index=threshold;index<=this->max;index++){
                        sigmaB += (index - mediaB)*(index - mediaB)*this->histogram[index]/(double)countB;
                }

                //3. Calcular la varianza inter clases

                double sigmaAB = wA*wB*(mediaA-mediaB)*(mediaA-mediaB);
                double criteria = (double)sigmaAB/(double)sigmaT;

                return criteria;

        }