16 double vtkOtsu::calculateOptimalThreshold(vtkImageData* volume, vtkImageData *region){
18 this->volume = volume;
19 this->region = region;
23 double valorEnRegion = 0.0;
31 this->volume->GetExtent(extent);
35 //1. Numero de voxeles
40 //2. Inicializar el histograma y criterios
41 this->histogram = new int[1000];
42 this->criterias = new double[1000];
43 for (index = 0; index <1000; index++){
44 this->histogram[index] = 0;
45 this->criterias[index] = 0.0;
49 //3. Encontrar el maximo, el promedio y calcular el histograma
51 for(i=extent[0];i<=extent[1];i++){
52 for(j=extent[2];j<=extent[3];j++){
53 for(k=extent[4];k<=extent[5];k++){
54 valorEnRegion = this->region->GetScalarComponentAsDouble(i,j,k,0);
55 if (valorEnRegion == 0.0){
56 valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
57 if (valor > this->max){
61 this->histogram[index]++;
69 media = (double)media/(double)this->voxelCount;
72 for(index=0;index<=max;index++){
73 this->sigmaT += (index - media)*(index - media)* this->histogram[index] / (double)this->voxelCount;
77 //calcular criterio para cada posible threshold e ir seleccionando el maximo
78 int optimalThreshold = 0;
79 double optimalCriteria = 0.0;
81 for (index=0;index<=max;index++){
82 criterias[index] = this->getCriteria(index);
83 if (criterias[index] > optimalCriteria){
84 optimalCriteria = criterias[index];
85 optimalThreshold = index;
88 return optimalThreshold;
94 double vtkOtsu::getCriteria(double threshold){
99 double valorEnRegion = 0.0;
101 this->volume->GetExtent(extent);
103 //1. Encontrar la frecuencia y la media de las clases A y B
112 for(i=extent[0];i<=extent[1];i++){
113 for(j=extent[2];j<=extent[3];j++){
114 for(k=extent[4];k<=extent[5];k++){
115 valorEnRegion = this->region->GetScalarComponentAsDouble(i,j,k,0);
116 if (valorEnRegion == 0.0){
117 valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
118 if (valor < threshold){
132 mediaA = (double)mediaA / (double)countA;
133 wA = (double)countA / (double)this->voxelCount;
135 mediaB = (double)mediaB / (double)countB;
136 wB = (double)countB / (double)this->voxelCount;
139 //2. Encontrar la varianza de las clases A y B
146 for(index=0;index<=threshold;index++){
147 sigmaA += (index - mediaA)*(index - mediaA)*this->histogram[index]/(double)countA;
151 for(index=threshold;index<=this->max;index++){
152 sigmaB += (index - mediaB)*(index - mediaB)*this->histogram[index]/(double)countB;
155 //3. Calcular la varianza inter clases
157 double sigmaAB = wA*wB*(mediaA-mediaB)*(mediaA-mediaB);
158 double criteria = (double)sigmaAB/(double)sigmaT;