/*# ---------------------------------------------------------------------
#
# Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
#                        pour la Sant�)
# Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
# Previous Authors : Laurent Guigues, Jean-Pierre Roux
# CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
#
#  This software is governed by the CeCILL-B license under French law and
#  abiding by the rules of distribution of free software. You can  use,
#  modify and/ or redistribute the software under the terms of the CeCILL-B
#  license as circulated by CEA, CNRS and INRIA at the following URL
#  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
#  or in the file LICENSE.txt.
#
#  As a counterpart to the access to the source code and  rights to copy,
#  modify and redistribute granted by the license, users are provided only
#  with a limited warranty  and the software's author,  the holder of the
#  economic rights,  and the successive licensors  have only  limited
#  liability.
#
#  The fact that you are presently reading this means that you have had
#  knowledge of the CeCILL-B license and that you accept its terms.
# ------------------------------------------------------------------------ */

#include "vtkOtsuImageData.h"

	 vtkOtsuImageData::vtkOtsuImageData(){
		 volume = NULL;
		 P = NULL;

	}

	 vtkOtsuImageData::~vtkOtsuImageData(){
		
	}

	int vtkOtsuImageData::calculateOptimalThreshold(vtkImageData* voi){
		
		int extent[6];
		int static TAM = 10000;
		
		
		this->volume = voi; //VOI
		this->volume->GetExtent(extent);

		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]);

					
		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
		
		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++){
					
										
					valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
					if (valor > max){
						max = (int)valor;
					}
					
					index  = (int)valor;
					this->H[index]++;
					voxelCount++;
					
				}
			}
		}

		printf("otsu: voxelCount = %d\n",voxelCount);
		
		//3. Normalizacion de la funcion P
		for(index=0;index<=max;index++){
			this->P[index] = (double)this->H[index] / (double)voxelCount;
			//printf("H[%d] = %d\n",index, H[index]);
		}




		//FILE *logger = NULL;
	
		//if ((logger = freopen("c:\\Diego\\Andes\\Tesis\\ITK\\FM_3D\\bin\\Debug\\thotsu.txt","w", stdout)) == NULL)
			//exit(-1);

		//printf("THRESHOLD; VARIANZA INTERCLASE\n");


		//4. Recorrer los diferentes niveles de intensidad

		for(T=0;T<=max;T++){

			//printf("T= %d, ",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;
			//printf("VI = %.2f\n",SI[T]);
		}
		
	    int optimalThreshold	= 0;
		double optimalCriteria	= VTK_LARGE_FLOAT;

	

		for (int index2=0;index2<=max;index2++){
			
			if (SI[index2] < optimalCriteria){
				optimalCriteria = SI[index2];
				optimalThreshold = index2;
				//printf("ENTRO %d %d\n",index2,optimalThreshold);
			}
			//printf("SI[%d] = %.2f, optimo = %d\n",index2,SI[index2],optimalThreshold);
		}
		//fclose(logger);
		//logger = NULL;

		printf(">>>%d<<<\n", optimalThreshold);
	
		return optimalThreshold;
	}