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

   1 /*# ---------------------------------------------------------------------
   2 #
   3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
   4 #                        pour la Sant�)
   5 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
   6 # Previous Authors : Laurent Guigues, Jean-Pierre Roux
   7 # CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
   8 #
   9 #  This software is governed by the CeCILL-B license under French law and
  10 #  abiding by the rules of distribution of free software. You can  use,
  11 #  modify and/ or redistribute the software under the terms of the CeCILL-B
  12 #  license as circulated by CEA, CNRS and INRIA at the following URL
  13 #  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  14 #  or in the file LICENSE.txt.
  15 #
  16 #  As a counterpart to the access to the source code and  rights to copy,
  17 #  modify and redistribute granted by the license, users are provided only
  18 #  with a limited warranty  and the software's author,  the holder of the
  19 #  economic rights,  and the successive licensors  have only  limited
  20 #  liability.
  21 #
  22 #  The fact that you are presently reading this means that you have had
  23 #  knowledge of the CeCILL-B license and that you accept its terms.
  24 # ------------------------------------------------------------------------ */
  25
  26 #include "vtkOtsuImageData.h"
  27
  28          vtkOtsuImageData::vtkOtsuImageData(){
  29                  volume = NULL;
  30                  P = NULL;
  31
  32         }
  33
  34          vtkOtsuImageData::~vtkOtsuImageData(){
  35
  36         }
  37
  38         int vtkOtsuImageData::calculateOptimalThreshold(vtkImageData* voi){
  39
  40                 int extent[6];
  41                 int static TAM = 10000;
  42
  43
  44                 this->volume = voi; //VOI
  45                 this->volume->GetExtent(extent);
  46
  47                 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]);
  48
  49
  50                 double valor = 0.0;
  51                 int index = 0;
  52                 int T = 0;
  53                 int max=0;
  54                 int voxelCount = 0;
  55
  56                 //1. Inicializacion
  57                 this->H = new int[TAM];
  58                 this->P = new double[TAM];
  59                 this->SI = new double[TAM];
  60                 for (index = 0; index <TAM; index++){
  61                         this->H[index] = 0;
  62                         this->P[index] = 0.0;   //probabilidad acumulada
  63                         this->SI[index] = 0.0;
  64                 }
  65
  66                 //2. Calcular el Histograma
  67
  68                 int i, j, k;
  69                 for(i=extent[0];i<=extent[1];i++){
  70                         for(j=extent[2];j<=extent[3];j++){
  71                                 for(k=extent[4];k<=extent[5];k++){
  72
  73
  74                                         valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
  75                                         if (valor > max){
  76                                                 max = (int)valor;
  77                                         }
  78
  79                                         index  = (int)valor;
  80                                         this->H[index]++;
  81                                         voxelCount++;
  82
  83                                 }
  84                         }
  85                 }
  86
  87                 printf("otsu: voxelCount = %d\n",voxelCount);
  88
  89                 //3. Normalizacion de la funcion P
  90                 for(index=0;index<=max;index++){
  91                         this->P[index] = (double)this->H[index] / (double)voxelCount;
  92                         //printf("H[%d] = %d\n",index, H[index]);
  93                 }
  94
  95
  96
  97
  98                 //FILE *logger = NULL;
  99
 100                 //if ((logger = freopen("c:\\Diego\\Andes\\Tesis\\ITK\\FM_3D\\bin\\Debug\\thotsu.txt","w", stdout)) == NULL)
 101                         //exit(-1);
 102
 103                 //printf("THRESHOLD; VARIANZA INTERCLASE\n");
 104
 105
 106                 //4. Recorrer los diferentes niveles de intensidad
 107
 108                 for(T=0;T<=max;T++){
 109
 110                         //printf("T= %d, ",T);
 111
 112                         //4.1 Calcular las probabilidades acumuladas de las clases A y B
 113                         double Q1 = 0.0;
 114                         double Q2 = 0.0;
 115                         for(index=0;index<=max;index++){
 116                                 if (index < T){
 117                                         Q1 = Q1 + this->P[index];
 118                                 }
 119                                 else{
 120                                         Q2 = Q2 + this->P[index];
 121                                 }
 122                         }
 123
 124                         //4.2 Calcular las medias de las clases A y B
 125                         double U1 = 0.0;
 126                         double U2 = 0.0;
 127                         for (index=0;index<=max;index++){
 128                                 if (index < T){
 129                                         U1 = U1 + ((double)(index * this->P[index])/(double)Q1);
 130                                 }
 131                                 else{
 132                                         U2 = U2 + ((double)(index * this->P[index])/(double)Q2);
 133                                 }
 134                         }
 135
 136
 137                         //4.3 Calcular las varianzas individuales
 138                         double S1 = 0.0;
 139                         double S2 = 0.0;
 140                         for (index=0;index<=max;index++){
 141                                 if (index < T){
 142                                         S1 = S1 + ((index - U1)*(index - U1)*(this->P[index]/(double)Q1));
 143                                 }
 144                                 else{
 145                                         S2 = S2 + ((index - U2)*(index - U2)*(this->P[index]/(double)Q2));
 146                                 }
 147                         }
 148
 149                         //4.4 Calcular SI: La varianza intraclases y seleccionar el valor minimo
 150                         this->SI[T] = Q1*S1 + Q2*S2;
 151                         //printf("VI = %.2f\n",SI[T]);
 152                 }
 153
 154             int optimalThreshold        = 0;
 155                 double optimalCriteria  = VTK_LARGE_FLOAT;
 156
 157
 158
 159                 for (int index2=0;index2<=max;index2++){
 160
 161                         if (SI[index2] < optimalCriteria){
 162                                 optimalCriteria = SI[index2];
 163                                 optimalThreshold = index2;
 164                                 //printf("ENTRO %d %d\n",index2,optimalThreshold);
 165                         }
 166                         //printf("SI[%d] = %.2f, optimo = %d\n",index2,SI[index2],optimalThreshold);
 167                 }
 168                 //fclose(logger);
 169                 //logger = NULL;
 170
 171                 printf(">>>%d<<<\n", optimalThreshold);
 172
 173                 return optimalThreshold;
 174         }
 175
 176
 177