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

   1 #include "vtkOtsu.h"
   2
   3          vtkOtsu::vtkOtsu(){
   4                  volume = NULL;
   5                  region = NULL;
   6
   7         }
   8
   9          vtkOtsu::~vtkOtsu(){
  10
  11         }
  12
  13
  14
  15
  16         double vtkOtsu::calculateOptimalThreshold(vtkImageData* volume, vtkImageData *region){
  17
  18                 this->volume = volume;
  19                 this->region = region;
  20
  21                 double valor = 0.0;
  22                 double media = 0.0;
  23                 double valorEnRegion = 0.0;
  24
  25
  26                 int extent[6];
  27
  28                 int index = 0;
  29
  30
  31                 this->volume->GetExtent(extent);
  32                 this->max=0.0;
  33                 this->sigmaT = 0.0;
  34
  35                 //1. Numero de voxeles
  36                 this->voxelCount = 0;
  37
  38
  39
  40                 //2. Inicializar el histograma y criterios
  41                 this->histogram = new int[1000];
  42                 this->criterias = new double[1000];
  43                 for (index = 0; index <1000; index++){
  44                         this->histogram[index] = 0;
  45                         this->criterias[index] = 0.0;
  46                 }
  47
  48
  49                 //3. Encontrar el maximo, el promedio y calcular el histograma
  50                 int i, j, k;
  51                 for(i=extent[0];i<=extent[1];i++){
  52                         for(j=extent[2];j<=extent[3];j++){
  53                                 for(k=extent[4];k<=extent[5];k++){
  54                                         valorEnRegion = this->region->GetScalarComponentAsDouble(i,j,k,0);
  55                                         if (valorEnRegion == 0.0){
  56                                                 valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
  57                                                 if (valor > this->max){
  58                                                         this->max = valor;
  59                                                 }
  60                                                 index  = (int)valor;
  61                                                 this->histogram[index]++;
  62                                                 media += valor;
  63                                                 this->voxelCount++;
  64                                         }
  65                                 }
  66                         }
  67                 }
  68
  69                 media = (double)media/(double)this->voxelCount;
  70
  71                 //calcula sigma total
  72                 for(index=0;index<=max;index++){
  73                         this->sigmaT += (index - media)*(index - media)* this->histogram[index] / (double)this->voxelCount;
  74                 }
  75
  76
  77                 //calcular criterio para cada posible threshold e ir seleccionando el maximo
  78             int optimalThreshold = 0;
  79                 double optimalCriteria = 0.0;
  80
  81                 for (index=0;index<=max;index++){
  82                         criterias[index] = this->getCriteria(index);
  83                         if (criterias[index] > optimalCriteria){
  84                                 optimalCriteria = criterias[index];
  85                                 optimalThreshold = index;
  86                         }
  87                 }
  88                 return optimalThreshold;
  89
  90
  91         }
  92
  93
  94         double vtkOtsu::getCriteria(double threshold){
  95
  96                 int extent[6];
  97                 int i,j,k;
  98                 double valor = 0.0;
  99                 double valorEnRegion = 0.0;
 100
 101                 this->volume->GetExtent(extent);
 102
 103                 //1. Encontrar la frecuencia y la media de las clases A y B
 104                 double mediaA=0.0;
 105                 double wA = 0.0;
 106                 int countA = 0;
 107
 108                 double mediaB=0.0;
 109                 double wB = 0.0;
 110                 int countB = 0;
 111
 112                 for(i=extent[0];i<=extent[1];i++){
 113                         for(j=extent[2];j<=extent[3];j++){
 114                                 for(k=extent[4];k<=extent[5];k++){
 115                                         valorEnRegion = this->region->GetScalarComponentAsDouble(i,j,k,0);
 116                                         if (valorEnRegion == 0.0){
 117                                                 valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
 118                                                 if (valor < threshold){
 119                                                         mediaA += valor;
 120                                                         countA++;
 121
 122                                                 }
 123                                                 else{
 124                                                         mediaB += valor;
 125                                                         countB++;
 126                                                 }
 127                                         }
 128                                 }
 129                         }
 130                 }
 131
 132                 mediaA = (double)mediaA / (double)countA;
 133                 wA = (double)countA / (double)this->voxelCount;
 134
 135                 mediaB = (double)mediaB / (double)countB;
 136                 wB = (double)countB / (double)this->voxelCount;
 137
 138
 139                 //2. Encontrar la varianza de las clases A y B
 140
 141                 double sigmaA=0.0;
 142                 double sigmaB=0.0;
 143
 144                 int index = 0;
 145
 146                 for(index=0;index<=threshold;index++){
 147                         sigmaA += (index - mediaA)*(index - mediaA)*this->histogram[index]/(double)countA;
 148                 }
 149
 150
 151                 for(index=threshold;index<=this->max;index++){
 152                         sigmaB += (index - mediaB)*(index - mediaB)*this->histogram[index]/(double)countB;
 153                 }
 154
 155                 //3. Calcular la varianza inter clases
 156
 157                 double sigmaAB = wA*wB*(mediaA-mediaB)*(mediaA-mediaB);
 158                 double criteria = (double)sigmaAB/(double)sigmaT;
 159
 160                 return criteria;
 161
 162         }