1 #include "vtkOtsuImageData.h"
3 vtkOtsuImageData::vtkOtsuImageData(){
9 vtkOtsuImageData::~vtkOtsuImageData(){
13 int vtkOtsuImageData::calculateOptimalThreshold(vtkImageData* voi){
16 int static TAM = 10000;
19 this->volume = voi; //VOI
20 this->volume->GetExtent(extent);
22 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]);
32 this->H = new int[TAM];
33 this->P = new double[TAM];
34 this->SI = new double[TAM];
35 for (index = 0; index <TAM; index++){
37 this->P[index] = 0.0; //probabilidad acumulada
38 this->SI[index] = 0.0;
41 //2. Calcular el Histograma
44 for(i=extent[0];i<=extent[1];i++){
45 for(j=extent[2];j<=extent[3];j++){
46 for(k=extent[4];k<=extent[5];k++){
49 valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
62 printf("otsu: voxelCount = %d\n",voxelCount);
64 //3. Normalizacion de la funcion P
65 for(index=0;index<=max;index++){
66 this->P[index] = (double)this->H[index] / (double)voxelCount;
67 //printf("H[%d] = %d\n",index, H[index]);
73 //FILE *logger = NULL;
75 //if ((logger = freopen("c:\\Diego\\Andes\\Tesis\\ITK\\FM_3D\\bin\\Debug\\thotsu.txt","w", stdout)) == NULL)
78 //printf("THRESHOLD; VARIANZA INTERCLASE\n");
81 //4. Recorrer los diferentes niveles de intensidad
85 //printf("T= %d, ",T);
87 //4.1 Calcular las probabilidades acumuladas de las clases A y B
90 for(index=0;index<=max;index++){
92 Q1 = Q1 + this->P[index];
95 Q2 = Q2 + this->P[index];
99 //4.2 Calcular las medias de las clases A y B
102 for (index=0;index<=max;index++){
104 U1 = U1 + ((double)(index * this->P[index])/(double)Q1);
107 U2 = U2 + ((double)(index * this->P[index])/(double)Q2);
112 //4.3 Calcular las varianzas individuales
115 for (index=0;index<=max;index++){
117 S1 = S1 + ((index - U1)*(index - U1)*(this->P[index]/(double)Q1));
120 S2 = S2 + ((index - U2)*(index - U2)*(this->P[index]/(double)Q2));
124 //4.4 Calcular SI: La varianza intraclases y seleccionar el valor minimo
125 this->SI[T] = Q1*S1 + Q2*S2;
126 //printf("VI = %.2f\n",SI[T]);
129 int optimalThreshold = 0;
130 double optimalCriteria = VTK_LARGE_FLOAT;
134 for (int index2=0;index2<=max;index2++){
136 if (SI[index2] < optimalCriteria){
137 optimalCriteria = SI[index2];
138 optimalThreshold = index2;
139 //printf("ENTRO %d %d\n",index2,optimalThreshold);
141 //printf("SI[%d] = %.2f, optimo = %d\n",index2,SI[index2],optimalThreshold);
146 printf(">>>%d<<<\n", optimalThreshold);
148 return optimalThreshold;