1 /*# ---------------------------------------------------------------------
3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
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
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.
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
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 # ------------------------------------------------------------------------ */
41 double vtkOtsu::calculateOptimalThreshold(vtkImageData* volume, vtkImageData *region){
43 this->volume = volume;
44 this->region = region;
48 double valorEnRegion = 0.0;
56 this->volume->GetExtent(extent);
60 //1. Numero de voxeles
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;
74 //3. Encontrar el maximo, el promedio y calcular el histograma
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){
86 this->histogram[index]++;
94 media = (double)media/(double)this->voxelCount;
97 for(index=0;index<=max;index++){
98 this->sigmaT += (index - media)*(index - media)* this->histogram[index] / (double)this->voxelCount;
102 //calcular criterio para cada posible threshold e ir seleccionando el maximo
103 int optimalThreshold = 0;
104 double optimalCriteria = 0.0;
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;
113 return optimalThreshold;
119 double vtkOtsu::getCriteria(double threshold){
124 double valorEnRegion = 0.0;
126 this->volume->GetExtent(extent);
128 //1. Encontrar la frecuencia y la media de las clases A y B
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){
157 mediaA = (double)mediaA / (double)countA;
158 wA = (double)countA / (double)this->voxelCount;
160 mediaB = (double)mediaB / (double)countB;
161 wB = (double)countB / (double)this->voxelCount;
164 //2. Encontrar la varianza de las clases A y B
171 for(index=0;index<=threshold;index++){
172 sigmaA += (index - mediaA)*(index - mediaA)*this->histogram[index]/(double)countA;
176 for(index=threshold;index<=this->max;index++){
177 sigmaB += (index - mediaB)*(index - mediaB)*this->histogram[index]/(double)countB;
180 //3. Calcular la varianza inter clases
182 double sigmaAB = wA*wB*(mediaA-mediaB)*(mediaA-mediaB);
183 double criteria = (double)sigmaAB/(double)sigmaT;