VMTK
equalize.cpp
1 /*
2  * equaliza.cpp: equaliza o histograma para aumentar a escala dinamica
3  * das intensidades
4  *
5  * Copyright (C) 2013 Wu Shin-Ting, FEEC, Unicamp
6  *
7  * This program is free software: you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation, either version 3 of the License, or
10  * any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program. If not, see <http://www.gnu.org/licenses/>.
19  *
20  */
21 #include <cmath>
22 #include <string.h>
23 #include "equalize.h"
24 
25 // Sugestão de leitura:
26 //http://en.wikipedia.org/wiki/Histogram_equalization
27 
28 void Equalize::EqualizeHistogram (int dimx, int dimy, int dimz,
29  unsigned short *volume, int nbitsalloc,
30  unsigned int umax, unsigned int **histogram)
31 {
32  // Computa o histograma do <volume>:
33  // frequencia de ocorrência de cada intensidade
34  int i;
35  int L = pow(2,nbitsalloc);
36  int *tmp = new int[L];
37  int *cdf = new int[L];
38 
39  // Initializa os vetores
40  memset(tmp, 0x00, L*sizeof(int));
41 
42  // Determina a frequencia de ocorrencia de cada intensidade
43  int intensidade, iz, iy, ix;
44  for (iz = 0; iz < dimz; iz++) {
45  for (iy =0; iy < dimy; iy++) {
46  for (ix =0; ix < dimx; ix++) {
47  intensidade = (int)(volume[iz*dimx*dimy+iy*dimx+ix]);
48  tmp[intensidade]++;
49  }
50  }
51  }
52 
53  int volSize = dimx*dimy*dimz;
54 
55  // Determina a funcao de distribuicao acumulada
56  // https://pt.wikipedia.org/wiki/Fun%C3%A7%C3%A3o_distribui%C3%A7%C3%A3o_acumulada
57  cdf[0] = tmp[0];
58  int cdfmin = -1;
59 
60  if (cdf[0]) cdfmin = cdf[0];
61  for(i = 1; i < L; i++) {
62  cdf[i] = cdf[i-1] + tmp[i];
63  if (cdfmin == -1 && cdf[i]) cdfmin = cdf[i];
64  }
65 
66  delete [] tmp;
67 
68  // Mapeamento de intensidades
69  *histogram = new unsigned int[umax];
70  memset(*histogram, 0x00, umax*sizeof(unsigned int));
71 
72  for (i=0; i<umax; i++) {
73  (*histogram)[i] = (unsigned int)(((1.0*(cdf[i]-cdfmin))/(volSize-cdfmin))*(L-1));
74  }
75  delete [] cdf;
76 }