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 # ------------------------------------------------------------------------ */
26 #include "vtkOtsuImageData.h"
28 vtkOtsuImageData::vtkOtsuImageData(){
34 vtkOtsuImageData::~vtkOtsuImageData(){
38 int vtkOtsuImageData::calculateOptimalThreshold(vtkImageData* voi){
41 int static TAM = 10000;
44 this->volume = voi; //VOI
45 this->volume->GetExtent(extent);
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]);
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++){
62 this->P[index] = 0.0; //probabilidad acumulada
63 this->SI[index] = 0.0;
66 //2. Calcular el Histograma
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++){
74 valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
87 printf("otsu: voxelCount = %d\n",voxelCount);
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]);
98 //FILE *logger = NULL;
100 //if ((logger = freopen("c:\\Diego\\Andes\\Tesis\\ITK\\FM_3D\\bin\\Debug\\thotsu.txt","w", stdout)) == NULL)
103 //printf("THRESHOLD; VARIANZA INTERCLASE\n");
106 //4. Recorrer los diferentes niveles de intensidad
110 //printf("T= %d, ",T);
112 //4.1 Calcular las probabilidades acumuladas de las clases A y B
115 for(index=0;index<=max;index++){
117 Q1 = Q1 + this->P[index];
120 Q2 = Q2 + this->P[index];
124 //4.2 Calcular las medias de las clases A y B
127 for (index=0;index<=max;index++){
129 U1 = U1 + ((double)(index * this->P[index])/(double)Q1);
132 U2 = U2 + ((double)(index * this->P[index])/(double)Q2);
137 //4.3 Calcular las varianzas individuales
140 for (index=0;index<=max;index++){
142 S1 = S1 + ((index - U1)*(index - U1)*(this->P[index]/(double)Q1));
145 S2 = S2 + ((index - U2)*(index - U2)*(this->P[index]/(double)Q2));
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]);
154 int optimalThreshold = 0;
155 double optimalCriteria = VTK_LARGE_FLOAT;
159 for (int index2=0;index2<=max;index2++){
161 if (SI[index2] < optimalCriteria){
162 optimalCriteria = SI[index2];
163 optimalThreshold = index2;
164 //printf("ENTRO %d %d\n",index2,optimalThreshold);
166 //printf("SI[%d] = %.2f, optimo = %d\n",index2,SI[index2],optimalThreshold);
171 printf(">>>%d<<<\n", optimalThreshold);
173 return optimalThreshold;