#include "vtkOtsuSphereSource.h"

	 vtkOtsuSphereSource::vtkOtsuSphereSource(){
		 volume = NULL;
		 region = NULL;
		 result = NULL;
		 P = NULL;

	}

	 vtkOtsuSphereSource::~vtkOtsuSphereSource(){
		
	}

	double vtkOtsuSphereSource::calculateOptimalThreshold(vtkImageData* voi){
		
		int extent[6];
		double centro[3];
		int static TAM = 2000;

			
		
		this->volume = voi; //VOI
		this->region = region; //JF
		this->volume->GetExtent(extent);

		centro[0] = (extent[0] + extent[1]) / 2;
		centro[1] = (extent[2] + extent[3]) / 2;
		centro[2] = (extent[4] + extent[5]) / 2;

		double radio = (extent[1] - extent[0]) /(double)2.0;

		result = vtkImageData::New();
     	result->DeepCopy(volume);
		
		double valorEnRegion = 0.0;
		double valor = 0.0;
		int index = 0;
		int T = 0;
		int max=0;
		int voxelCount = 0;
	
		//1. Inicializacion
		this->H = new int[TAM];
		this->P = new double[TAM];
		this->SI = new double[TAM];
		for (index = 0; index <TAM; index++){
			this->H[index] = 0;
			this->P[index] = 0.0;   //probabilidad acumulada
			this->SI[index] = 0.0;
		}
		
		//2. Calcular el Histograma
		double *ptr;
		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++){
					
					ptr = (double *)this->result->GetScalarPointer(i,j,k);

					// verificar si esta dentro de la esfera
					if(  ((i-centro[0])*(i-centro[0]) + 
						 (j-centro[1])*(j-centro[1]) +
						 (k-centro[2])*(k-centro[2])) <= (radio*radio)){

						
							valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
							if (valor > max){
								max = (int)valor;
							}
							index  = (int)valor;
							this->H[index]++;
							voxelCount++;
						
						
					}
					else {
						*ptr=0.0;
					}
					
				}
			}
		}
		
		//3. Normalizacion de la funcion P
		for(index=0;index<=max;index++){
			this->P[index] = (double)this->H[index] / (double)voxelCount;
		}


		//4. Recorrer los diferentes niveles de intensidad

		for(T=0;T<=max;T++){
			
			//4.1 Calcular las probabilidades acumuladas de las clases A y B
			double Q1 = 0.0;
			double Q2 = 0.0;
			for(index=0;index<=max;index++){
				if (index < T){
					Q1 = Q1 + this->P[index];
				}
				else{
					Q2 = Q2 + this->P[index];
				}
			}

			//4.2 Calcular las medias de las clases A y B
			double U1 = 0.0;
			double U2 = 0.0;
			for (index=0;index<=max;index++){
				if (index < T){
					U1 = U1 + ((double)(index * this->P[index])/(double)Q1);
				}
				else{
					U2 = U2 + ((double)(index * this->P[index])/(double)Q2);
				}
			}

			
			//4.3 Calcular las varianzas individuales
			double S1 = 0.0;
			double S2 = 0.0;
			for (index=0;index<=max;index++){
				if (index < T){
					S1 = S1 + ((index - U1)*(index - U1)*(this->P[index]/(double)Q1));
				}
				else{
					S2 = S2 + ((index - U2)*(index - U2)*(this->P[index]/(double)Q2));
				}
			}

			//4.4 Calcular SI: La varianza intraclases y seleccionar el valor minimo
			this->SI[T] = Q1*S1 + Q2*S2;
		}
		
	    int optimalThreshold	= 0;
		double optimalCriteria	= VTK_LARGE_FLOAT;

		for (index=0;index<=max;index++){
			if (SI[index] < optimalCriteria){
				optimalCriteria = SI[index];
				optimalThreshold = index;
			}
		}
		return optimalThreshold;
	}




	vtkImageData* vtkOtsuSphereSource::getImageForSegmentation(){
		return this->result;
	}