#include "vtkOtsu.h"

	 vtkOtsu::vtkOtsu(){
		 volume = NULL;
		 region = NULL;
	
	}

	 vtkOtsu::~vtkOtsu(){
		
	}




	double vtkOtsu::calculateOptimalThreshold(vtkImageData* volume, vtkImageData *region){

		this->volume = volume;
		this->region = region;
	
		double valor = 0.0;
		double media = 0.0;
		double valorEnRegion = 0.0;
	    
	
		int extent[6];
		
		int index = 0;
		
		
		this->volume->GetExtent(extent);
		this->max=0.0;
		this->sigmaT = 0.0;
		
		//1. Numero de voxeles
		this->voxelCount = 0;
		


		//2. Inicializar el histograma y criterios
		this->histogram = new int[1000];
		this->criterias = new double[1000];
		for (index = 0; index <1000; index++){
			this->histogram[index] = 0;
			this->criterias[index] = 0.0;
		}
		
		
		//3. Encontrar el maximo, el promedio y calcular el histograma
		int i, j, k;
		for(i=extent[0];i<=extent[1];i++){
			for(j=extent[2];j<=extent[3];j++){
				for(k=extent[4];k<=extent[5];k++){
					valorEnRegion = this->region->GetScalarComponentAsDouble(i,j,k,0);
					if (valorEnRegion == 0.0){
						valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
						if (valor > this->max){
							this->max = valor;
						}
						index  = (int)valor;
						this->histogram[index]++;
						media += valor;
						this->voxelCount++;
					}
				}
			}
		}
		
		media = (double)media/(double)this->voxelCount;

		//calcula sigma total
		for(index=0;index<=max;index++){
			this->sigmaT += (index - media)*(index - media)* this->histogram[index] / (double)this->voxelCount;
		}


		//calcular criterio para cada posible threshold e ir seleccionando el maximo
	    int optimalThreshold = 0;
		double optimalCriteria = 0.0;

		for (index=0;index<=max;index++){
			criterias[index] = this->getCriteria(index);
			if (criterias[index] > optimalCriteria){
				optimalCriteria = criterias[index];
				optimalThreshold = index;
			}
		}
		return optimalThreshold;

		
	}


	double vtkOtsu::getCriteria(double threshold){
		
		int extent[6];
		int i,j,k;
		double valor = 0.0;
		double valorEnRegion = 0.0;

		this->volume->GetExtent(extent);
		
		//1. Encontrar la frecuencia y la media de las clases A y B
		double mediaA=0.0;
		double wA = 0.0;
		int countA = 0;
		
		double mediaB=0.0;
		double wB = 0.0;
		int countB = 0;

		for(i=extent[0];i<=extent[1];i++){
			for(j=extent[2];j<=extent[3];j++){
				for(k=extent[4];k<=extent[5];k++){
					valorEnRegion = this->region->GetScalarComponentAsDouble(i,j,k,0);
					if (valorEnRegion == 0.0){
						valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
						if (valor < threshold){
							mediaA += valor;
							countA++;

						}
						else{
							mediaB += valor;
							countB++;
						}
					}
				}
			}
		}

		mediaA = (double)mediaA / (double)countA;
		wA = (double)countA / (double)this->voxelCount;

		mediaB = (double)mediaB / (double)countB;
		wB = (double)countB / (double)this->voxelCount;


		//2. Encontrar la varianza de las clases A y B

		double sigmaA=0.0;
		double sigmaB=0.0;

		int index = 0;
		
		for(index=0;index<=threshold;index++){
			sigmaA += (index - mediaA)*(index - mediaA)*this->histogram[index]/(double)countA;
		}

		
		for(index=threshold;index<=this->max;index++){
			sigmaB += (index - mediaB)*(index - mediaB)*this->histogram[index]/(double)countB;
		}

		//3. Calcular la varianza inter clases

		double sigmaAB = wA*wB*(mediaA-mediaB)*(mediaA-mediaB);
		double criteria = (double)sigmaAB/(double)sigmaT;

		return criteria;

	}