www.pudn.com > source.zip > gabor.cpp
#include "gui.h"
void ImageOperations::Gabor(){
/* Set default values for options
*/
wxString txt;
int counter = 0;
num_rows = width;
num_cols = height;
do_stat = FALSE;
save_filter = FALSE;
no_pgm = FALSE;
save_freq_resp = FALSE;
save_resp = FALSE;
select_filters = FALSE;
nogauss_lowpass = FALSE;
nononlin_tanh = FALSE;
alpha = DEFAULT_ALPHA;
LPsigmafact = DEFAULT_SIGMA_FACTOR;
satisfy_COD = DEFAULT_COD;
/* Scale input image to [-1,1] and set the mean value to zero in order
* to avoid filter response to the mean gray value.
*/
image = gray_to_real(GrayImage, NULL, num_rows, num_cols, num_cols);
/* Mirror pad image in order to evitate aliasing in convolution.
* mirror_pad_real() allocates a fftw_real type array for the padded
* input image which will also hold its Fourier-transform. Since the
* transform has fftw_complex type, we have to allocate a bigger
* array, namely p_num_rows X p_ft_num_cols. For more deatails, see
* the FFTW documentation (http://www.fftw.org) under Sections 2.4
* (Real Multi-dimensional Transforms Tutorial) and 3.5 (Real
* Multi-dimensional Transforms Reference).
*/
row_pad = num_rows / 2;
col_pad = num_cols / 2;
p_num_rows = num_rows + 2 * row_pad;
p_num_cols = num_cols + 2 * col_pad;
p_ft_num_cols = 2 * (p_num_cols/2 + 1);
p_image = mirror_pad_real(image, NULL, num_rows, num_cols, num_cols,
p_ft_num_cols, row_pad, col_pad);
/* alias for accessing complex transform output.
*/
p_ft_image = (fftw_complex*) p_image;
/* Allocate filter of the same size as p_image. response is an alias to
* access the response image after inverse FFT.
*/
filter = (fftw_complex *) Alloc_dvector(p_num_rows * p_ft_num_cols);
response = (fftw_real *) filter;
/* Compute Fourier transform of the image in-place.
* We use the FFTW library available at http://www.fftw.org .
* Details about the functions can be found in the documentation
* under Sections 2.4 (Real Multi-dimensional Transforms Tutorial) and
* 3.5 (Real Multi-dimensional Transforms Reference).
*/
p = rfftw2d_create_plan(p_num_rows, p_num_cols, FFTW_REAL_TO_COMPLEX,
FFTW_ESTIMATE | FFTW_IN_PLACE);
pinv = rfftw2d_create_plan(p_num_rows, p_num_cols, FFTW_COMPLEX_TO_REAL,
FFTW_ESTIMATE | FFTW_IN_PLACE);
rfftwnd_one_real_to_complex(p, p_image, NULL);
/* Scale factor for convolution. FFTW transforms are unnormalized, we
* now normalize it.
*/
scale = 1.0/((double)p_num_rows * p_num_cols);
for (i = 0; i < p_num_rows * p_ft_num_cols/2; i++) {
p_ft_image[i].re *= scale;
p_ft_image[i].im *= scale;
}
/* Number of radial frequencies and num. of lowest radial frequencies to skip.
*/
/* modify 1 */
num_rad_freq_double = (log((double)num_cols / 2 )/log(2));
num_rad_freq = (int) num_rad_freq_double;
skip_low = (num_rad_freq - NUM_RAD_FREQ_SKIP > 0?
NUM_RAD_FREQ_SKIP : 0);
num_rad_freq -= skip_low;
E = Alloc_dvector(num_rad_freq * NUM_ORIENT);
Esorted = Alloc_ivector(num_rad_freq * NUM_ORIENT);
U = Alloc_dvector(num_rad_freq);
SigmaU = Alloc_dvector(num_rad_freq);
SigmaV = Alloc_dvector(num_rad_freq);
/* Normalizing factor for U
*/
U_norm_fact = ((double)num_cols - 1.0) / SQ((double)num_cols);
/* Set up Gabor filter parameters U, V and Theta.
*/
for (i = 0; i < num_rad_freq; i++) {
U[i] = pow(2.0, (double)(i + skip_low)) * M_SQRT2 * U_norm_fact;
SigmaU[i] = DEFAULT_SIGMA_U_FACTOR * U[i];
SigmaV[i] = DEFAULT_SIGMA_V_FACTOR * U[i];
}
for (i = 0; i < NUM_ORIENT; i++) Theta[i] = (double)i * M_PI_4;
/* Compute one response-image per filter
*/
Esum = 0.0;
realMatrix = new fftw_real**[num_rad_freq * NUM_ORIENT];
for (i=0; i<(num_rad_freq * NUM_ORIENT); ++i){
realMatrix[i] = new fftw_real*[width];
for (j = 0; j < width; ++j){
realMatrix[i][j] = new fftw_real[height];
}
}
for (i = 0; i < num_rad_freq; i++){
for (j = 0; j < NUM_ORIENT; j++) {
ft_even_gabor(filter, p_num_rows, p_num_cols, U[i], SigmaU[i], SigmaV[i],
Theta[j]);
if (save_filter) {
gray **tmp;
tmp = halfcomplex_to_gray(filter, NULL, p_num_rows, p_num_cols,
p_ft_num_cols/2);
sprintf(fname, "gabor.%uX%u.%u.%u.bmp\0", num_cols, num_rows, i, j);
tmp_image = new wxImage(width, height, dim(tmp));
tmp_image->SaveFile(fname, wxBITMAP_TYPE_BMP);
}
/* Filter convolution. It is computed in-place in filter.
*/
l = i*NUM_ORIENT+j;
E[l] = 0.0;
for (row = 0; row < p_num_rows; row++)
for (col = 0; col < p_ft_num_cols/2; col++) {
k = row * p_ft_num_cols/2 + col;
filter[k].im = (p_ft_image[k].im * filter[k].re);
filter[k].re = (p_ft_image[k].re * filter[k].re);
if (row >= row_pad && row < num_rows + row_pad &&
col >= col_pad && col < num_cols + col_pad)
E[l] += (SQ(filter[k].re) + SQ(filter[k].im));
}
Esum += E[l];
if (save_freq_resp) {
/* Save frequency domain filter response.
*/
sprintf(fname, "%s.%u.%u.Gabor1\0", "Kimenet", i, j);
save_halfcomplexmatrix(fname, p_num_rows, p_ft_num_cols/2, 0, filter);
}
/* Invers FFT in-place in filter.
*/
rfftwnd_one_complex_to_real(pinv, filter, NULL);
/* Extract filtered image into 'image' and scale to [-1,1]
*/
image = extract_realsubimage(response, image, num_rows, num_cols,
p_ft_num_cols, num_cols, row_pad, col_pad);
image = real_scaling(image, image, num_rows, num_cols,
num_cols, num_cols);
/* modify 2 */
/* Save response-image, load them later if needed when doing feature
extraction*/
sprintf(fname, "%s.%u.%u.gabor2\0", ".\\tmp\\Kimenet", i, j);
for (int i2 = 0; i2 < num_rows; i2++){
for (int j2 = 0; j2 < num_cols; j2++){
realMatrix[i*num_rad_freq+j][i2][j2] = image[i2*(num_cols)+j2];
}
}
// save_realmatrix(fname, num_rows, num_cols, 0, image);
if (no_pgm) {
/* Scale and save filtered image in PGM format
*/
in = NULL;
in = real_to_gray(image, in, num_rows, num_cols, num_cols);
sprintf(fname, "%s.%u.%u.gabor.bmp\0", "Kimenet", i, j);
//save_image(fname, NULL, 1, num_cols, num_rows, &in);
tmp_image = new wxImage(num_rows, num_cols , dim(in));
tmp_image->SaveFile(fname, wxBITMAP_TYPE_BMP);
}
}
}//end of loop//
/* filter selection */
/* Compute approximate COD (coefficient of determination) for each filtered
* image.
*/
index = (int *)malloc(200*sizeof(int));
for (i=0; i < 100; i++){index[i] = -1;}
((MyFrame *)frame)->gaussians->AppendText(txt);
for (i = 0; i < num_rad_freq * NUM_ORIENT; i++) {
Esorted[i] = i;
E[i] /= Esum;
}
/* Sort COD's in descending order and print them on stdout.
*/
qsort(Esorted, num_rad_freq * NUM_ORIENT, sizeof(int), Ecmp);
Esum = 0.0;
last_filter = 0;
printf("Approximate COD (coefficient of determination) for each filtered "\
"image in descending order:\n");
for (i = num_rad_freq * NUM_ORIENT-1; i >=0 ; i--) {
// printf("%u %u E[%u][%u] = %e (Total = %e)\n", Esorted[i]/NUM_ORIENT,Esorted[i]%NUM_ORIENT,Esorted[i]/NUM_ORIENT,
// Esorted[i]%NUM_ORIENT, E[Esorted[i]], Esum+=E[Esorted[i]]);
Esum+=E[Esorted[i]];
// txt = "";
// txt << Esorted[i]/NUM_ORIENT<< " "<< Esorted[i]%NUM_ORIENT <<" E["
// << Esorted[i]/NUM_ORIENT<< "]["<< Esorted[i]%NUM_ORIENT << "] = "
// << E[Esorted[i]] << " (Total = " << Esum;
// ((MyFrame *)frame)->gaussians->AppendText(txt);
// ((MyFrame *)frame)->gaussians->AppendText("\n");
if (Esum >= satisfy_COD) {
last_filter = i;
break;
}
}
///
//last_filter = 0;
/* Feature extraction for the selected filtered image */
Texture = new gray**[num_rad_freq * NUM_ORIENT];//-last_filter];
for (i=0; i<(num_rad_freq * NUM_ORIENT); ++i){//-last_filter];
Texture[i] = new gray*[width];
for (j = 0; j < width; ++j){
Texture[i][j] = new gray[height];
}
}
no_g_feature = num_rad_freq * NUM_ORIENT-last_filter;
((MyFrame *)frame)->Set_Gabor(num_rad_freq * NUM_ORIENT);
int cntr = 0;
//for (f = num_rad_freq * NUM_ORIENT-1; f>=last_filter ; f--) {
for (f = num_rad_freq * NUM_ORIENT-1; f>=0 ; f--) {
i = Esorted[f]/NUM_ORIENT;
j = Esorted[f]%NUM_ORIENT;
index[2*counter] = i; index[2*counter+1] = j;
counter++;
sprintf(fname, "%s.%u.%u.gabor2\0", ".\\tmp\\Kimenet", i, j);
for (int i2 = 0; i2 < num_rows; i2++){
for (int j2 = 0; j2 < num_cols; j2++){
image[i2*(num_cols)+j2] = realMatrix[i*num_rad_freq+j][i2][j2];
}
}
// image = load_realmatrix(fname, &num_rows, &num_cols, 0);
if(!nononlin_tanh) {
/* Nonlinear transform of filtered image. We use herein a
* bounded nonlinearity similar to the sigmoidal activation
* function used in artificial neural networks. Original
* equations can be found in:
* A. K. Jain and F. Farrokhnia:Unsupervised texture Segmentation
* Using Gabor Filters. Pattern Recognition Vol. 24, No. 12,
* pp 1167-1186, 1991.
*/
for (row = 0; row < num_rows; row++)
for (col = 0; col < num_cols; col++)
image[row*num_cols + col] = ABS(tanh(alpha * image[row*num_cols + col]));
}
if (!nogauss_lowpass) {
/* Apply a Gaussian lowpass filter.
* Filtering is done in the Fourier domain.
*/
double u, v, gauss, sigma, halfInvSigmaSq;
sigma = (U[i] == 0.0? LPsigmafact : LPsigmafact / U[i]);
printf("U[%u]=%e, sigma=%e\n", i, U[i], sigma);
/* Mirror pad image in 'response'.
*/
response = mirror_pad_real(image, response, num_rows, num_cols,
num_cols, p_ft_num_cols, row_pad, col_pad);
/* FFT padded image.
*/
rfftwnd_one_real_to_complex(p, response, NULL);
/* Convert sigma from spatial to frequency domain and
* compute convolution.
*/
sigma = 1.0 / (2.0 * M_PI * sigma);
halfInvSigmaSq = 0.5 / SQ(sigma);
for (row = 0; row < p_num_rows; row++) {
u = (row < p_num_rows/2? (double)row/p_num_rows :
(double) -(p_num_rows-row)/p_num_rows);
for (col = 0; col < p_ft_num_cols/2; col++) {
v = (double)col/p_num_cols;
gauss = exp(-(SQ(u) + SQ(v)) * halfInvSigmaSq);
k = row * p_ft_num_cols/2 + col;
filter[k].re *= gauss;
filter[k].im *= gauss;
}
}
/* IFFT response and extract filtered image.
*/
rfftwnd_one_complex_to_real(pinv, filter, NULL);
image = extract_realsubimage(response, image, num_rows, num_cols,
p_ft_num_cols, num_cols, row_pad, col_pad);
}
if (save_feature_pgm) {
/* Scale and save feature image in PGM format
*/
in = NULL;
in = real_to_gray(image, in, num_rows, num_cols, num_cols);
sprintf(fname, "%s.%u.%u.gabor.feature.bmp\0", ".\\tmp\\Kimenet", i, j);
//save_image(fname, NULL, 1, num_cols, num_rows, &in);
tmp_image = new wxImage(num_rows, num_cols, dim(in));
//tmp_image->SaveFile(fname, wxBITMAP_TYPE_BMP);
unsigned char *tmp_data = tmp_image->GetData();
for (i=0; i