1 #include "vtkOtsuSphereSource.h"
3 vtkOtsuSphereSource::vtkOtsuSphereSource(){
11 vtkOtsuSphereSource::~vtkOtsuSphereSource(){
15 double vtkOtsuSphereSource::calculateOptimalThreshold(vtkImageData* voi){
19 int static TAM = 2000;
23 this->volume = voi; //VOI
24 this->region = region; //JF
25 this->volume->GetExtent(extent);
27 centro[0] = (extent[0] + extent[1]) / 2;
28 centro[1] = (extent[2] + extent[3]) / 2;
29 centro[2] = (extent[4] + extent[5]) / 2;
31 double radio = (extent[1] - extent[0]) /(double)2.0;
33 result = vtkImageData::New();
34 result->DeepCopy(volume);
36 double valorEnRegion = 0.0;
44 this->H = new int[TAM];
45 this->P = new double[TAM];
46 this->SI = new double[TAM];
47 for (index = 0; index <TAM; index++){
49 this->P[index] = 0.0; //probabilidad acumulada
50 this->SI[index] = 0.0;
53 //2. Calcular el Histograma
56 for(i=extent[0];i<=extent[1];i++){
57 for(j=extent[2];j<=extent[3];j++){
58 for(k=extent[4];k<=extent[5];k++){
60 ptr = (double *)this->result->GetScalarPointer(i,j,k);
62 // verificar si esta dentro de la esfera
63 if( ((i-centro[0])*(i-centro[0]) +
64 (j-centro[1])*(j-centro[1]) +
65 (k-centro[2])*(k-centro[2])) <= (radio*radio)){
68 valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
86 //3. Normalizacion de la funcion P
87 for(index=0;index<=max;index++){
88 this->P[index] = (double)this->H[index] / (double)voxelCount;
92 //4. Recorrer los diferentes niveles de intensidad
96 //4.1 Calcular las probabilidades acumuladas de las clases A y B
99 for(index=0;index<=max;index++){
101 Q1 = Q1 + this->P[index];
104 Q2 = Q2 + this->P[index];
108 //4.2 Calcular las medias de las clases A y B
111 for (index=0;index<=max;index++){
113 U1 = U1 + ((double)(index * this->P[index])/(double)Q1);
116 U2 = U2 + ((double)(index * this->P[index])/(double)Q2);
121 //4.3 Calcular las varianzas individuales
124 for (index=0;index<=max;index++){
126 S1 = S1 + ((index - U1)*(index - U1)*(this->P[index]/(double)Q1));
129 S2 = S2 + ((index - U2)*(index - U2)*(this->P[index]/(double)Q2));
133 //4.4 Calcular SI: La varianza intraclases y seleccionar el valor minimo
134 this->SI[T] = Q1*S1 + Q2*S2;
137 int optimalThreshold = 0;
138 double optimalCriteria = VTK_LARGE_FLOAT;
140 for (index=0;index<=max;index++){
141 if (SI[index] < optimalCriteria){
142 optimalCriteria = SI[index];
143 optimalThreshold = index;
146 return optimalThreshold;
152 vtkImageData* vtkOtsuSphereSource::getImageForSegmentation(){