lib/maracasVisuLib/src/kernel/vtkOtsuSphereSource.cxx

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