1 /*# ---------------------------------------------------------------------
3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
5 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
6 # Previous Authors : Laurent Guigues, Jean-Pierre Roux
7 # CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
9 # This software is governed by the CeCILL-B license under French law and
10 # abiding by the rules of distribution of free software. You can use,
11 # modify and/ or redistribute the software under the terms of the CeCILL-B
12 # license as circulated by CEA, CNRS and INRIA at the following URL
13 # http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
14 # or in the file LICENSE.txt.
16 # As a counterpart to the access to the source code and rights to copy,
17 # modify and redistribute granted by the license, users are provided only
18 # with a limited warranty and the software's author, the holder of the
19 # economic rights, and the successive licensors have only limited
22 # The fact that you are presently reading this means that you have had
23 # knowledge of the CeCILL-B license and that you accept its terms.
24 # ------------------------------------------------------------------------ */
26 #include "vtkOtsuSphereSource.h"
28 vtkOtsuSphereSource::vtkOtsuSphereSource(){
36 vtkOtsuSphereSource::~vtkOtsuSphereSource(){
40 double vtkOtsuSphereSource::calculateOptimalThreshold(vtkImageData* voi){
44 int static TAM = 2000;
48 this->volume = voi; //VOI
49 this->region = region; //JF
50 this->volume->GetExtent(extent);
52 centro[0] = (extent[0] + extent[1]) / 2;
53 centro[1] = (extent[2] + extent[3]) / 2;
54 centro[2] = (extent[4] + extent[5]) / 2;
56 double radio = (extent[1] - extent[0]) /(double)2.0;
58 result = vtkImageData::New();
59 result->DeepCopy(volume);
61 double valorEnRegion = 0.0;
69 this->H = new int[TAM];
70 this->P = new double[TAM];
71 this->SI = new double[TAM];
72 for (index = 0; index <TAM; index++){
74 this->P[index] = 0.0; //probabilidad acumulada
75 this->SI[index] = 0.0;
78 //2. Calcular el Histograma
81 for(i=extent[0];i<=extent[1];i++){
82 for(j=extent[2];j<=extent[3];j++){
83 for(k=extent[4];k<=extent[5];k++){
85 ptr = (double *)this->result->GetScalarPointer(i,j,k);
87 // verificar si esta dentro de la esfera
88 if( ((i-centro[0])*(i-centro[0]) +
89 (j-centro[1])*(j-centro[1]) +
90 (k-centro[2])*(k-centro[2])) <= (radio*radio)){
93 valor = this->volume->GetScalarComponentAsDouble(i,j,k,0);
111 //3. Normalizacion de la funcion P
112 for(index=0;index<=max;index++){
113 this->P[index] = (double)this->H[index] / (double)voxelCount;
117 //4. Recorrer los diferentes niveles de intensidad
121 //4.1 Calcular las probabilidades acumuladas de las clases A y B
124 for(index=0;index<=max;index++){
126 Q1 = Q1 + this->P[index];
129 Q2 = Q2 + this->P[index];
133 //4.2 Calcular las medias de las clases A y B
136 for (index=0;index<=max;index++){
138 U1 = U1 + ((double)(index * this->P[index])/(double)Q1);
141 U2 = U2 + ((double)(index * this->P[index])/(double)Q2);
146 //4.3 Calcular las varianzas individuales
149 for (index=0;index<=max;index++){
151 S1 = S1 + ((index - U1)*(index - U1)*(this->P[index]/(double)Q1));
154 S2 = S2 + ((index - U2)*(index - U2)*(this->P[index]/(double)Q2));
158 //4.4 Calcular SI: La varianza intraclases y seleccionar el valor minimo
159 this->SI[T] = Q1*S1 + Q2*S2;
162 int optimalThreshold = 0;
163 double optimalCriteria = VTK_LARGE_FLOAT;
165 for (index=0;index<=max;index++){
166 if (SI[index] < optimalCriteria){
167 optimalCriteria = SI[index];
168 optimalThreshold = index;
171 return optimalThreshold;
177 vtkImageData* vtkOtsuSphereSource::getImageForSegmentation(){