lib/maracasVisuLib/src/kernel/include/vtkOtsuImageData.cxx

   1 #include "vtkOtsuImageData.h"
   2
   3          vtkOtsuImageData::vtkOtsuImageData(){
   4                  volume = NULL;
   5                  P = NULL;
   6
   7         }
   8
   9          vtkOtsuImageData::~vtkOtsuImageData(){
  10
  11         }
  12
  13         int vtkOtsuImageData::calculateOptimalThreshold(vtkImageData* voi){
  14
  15                 int extent[6];
  16                 int static TAM = 10000;
  17
  18
  19                 this->volume = voi; //VOI
  20                 this->volume->GetExtent(extent);
  21
  22                 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]);
  23
  24
  25                 double valor = 0.0;
  26                 int index = 0;
  27                 int T = 0;
  28                 int max=0;
  29                 int voxelCount = 0;
  30
  31                 //1. Inicializacion
  32                 this->H = new int[TAM];
  33                 this->P = new double[TAM];
  34                 this->SI = new double[TAM];
  35                 for (index = 0; index <TAM; index++){
  36                         this->H[index] = 0;
  37                         this->P[index] = 0.0;   //probabilidad acumulada
  38                         this->SI[index] = 0.0;
  39                 }
  40
  41                 //2. Calcular el Histograma
  42
  43                 int i, j, k;
  44                 for(i=extent[0];i<=extent[1];i++){
  45                         for(j=extent[2];j<=extent[3];j++){
  46                                 for(k=extent[4];k<=extent[5];k++){
  47
  48
  49                                         valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
  50                                         if (valor > max){
  51                                                 max = (int)valor;
  52                                         }
  53
  54                                         index  = (int)valor;
  55                                         this->H[index]++;
  56                                         voxelCount++;
  57
  58                                 }
  59                         }
  60                 }
  61
  62                 printf("otsu: voxelCount = %d\n",voxelCount);
  63
  64                 //3. Normalizacion de la funcion P
  65                 for(index=0;index<=max;index++){
  66                         this->P[index] = (double)this->H[index] / (double)voxelCount;
  67                         //printf("H[%d] = %d\n",index, H[index]);
  68                 }
  69
  70
  71
  72
  73                 //FILE *logger = NULL;
  74
  75                 //if ((logger = freopen("c:\\Diego\\Andes\\Tesis\\ITK\\FM_3D\\bin\\Debug\\thotsu.txt","w", stdout)) == NULL)
  76                         //exit(-1);
  77
  78                 //printf("THRESHOLD; VARIANZA INTERCLASE\n");
  79
  80
  81                 //4. Recorrer los diferentes niveles de intensidad
  82
  83                 for(T=0;T<=max;T++){
  84
  85                         //printf("T= %d, ",T);
  86
  87                         //4.1 Calcular las probabilidades acumuladas de las clases A y B
  88                         double Q1 = 0.0;
  89                         double Q2 = 0.0;
  90                         for(index=0;index<=max;index++){
  91                                 if (index < T){
  92                                         Q1 = Q1 + this->P[index];
  93                                 }
  94                                 else{
  95                                         Q2 = Q2 + this->P[index];
  96                                 }
  97                         }
  98
  99                         //4.2 Calcular las medias de las clases A y B
 100                         double U1 = 0.0;
 101                         double U2 = 0.0;
 102                         for (index=0;index<=max;index++){
 103                                 if (index < T){
 104                                         U1 = U1 + ((double)(index * this->P[index])/(double)Q1);
 105                                 }
 106                                 else{
 107                                         U2 = U2 + ((double)(index * this->P[index])/(double)Q2);
 108                                 }
 109                         }
 110
 111
 112                         //4.3 Calcular las varianzas individuales
 113                         double S1 = 0.0;
 114                         double S2 = 0.0;
 115                         for (index=0;index<=max;index++){
 116                                 if (index < T){
 117                                         S1 = S1 + ((index - U1)*(index - U1)*(this->P[index]/(double)Q1));
 118                                 }
 119                                 else{
 120                                         S2 = S2 + ((index - U2)*(index - U2)*(this->P[index]/(double)Q2));
 121                                 }
 122                         }
 123
 124                         //4.4 Calcular SI: La varianza intraclases y seleccionar el valor minimo
 125                         this->SI[T] = Q1*S1 + Q2*S2;
 126                         //printf("VI = %.2f\n",SI[T]);
 127                 }
 128
 129             int optimalThreshold        = 0;
 130                 double optimalCriteria  = VTK_LARGE_FLOAT;
 131
 132
 133
 134                 for (int index2=0;index2<=max;index2++){
 135
 136                         if (SI[index2] < optimalCriteria){
 137                                 optimalCriteria = SI[index2];
 138                                 optimalThreshold = index2;
 139                                 //printf("ENTRO %d %d\n",index2,optimalThreshold);
 140                         }
 141                         //printf("SI[%d] = %.2f, optimo = %d\n",index2,SI[index2],optimalThreshold);
 142                 }
 143                 //fclose(logger);
 144                 //logger = NULL;
 145
 146                 printf(">>>%d<<<\n", optimalThreshold);
 147
 148                 return optimalThreshold;
 149         }
 150
 151
 152