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