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