www.pudn.com > source.zip > gaborutil.cpp
#include "gaborutil.h"
char ** pm_allocarray(int cols, int rows, int size)
{
char **its;
int i;
its = (char **) malloc(rows * sizeof(char *));
if (its == (char **) 0) {
fprintf(stderr, "%s: out of memory allocating an array\n", "demo");
return (char **) 0;
}
its[0] = (char *) malloc(rows * cols * size);
if (its[0] == (char *) 0) {
fprintf(stderr, "%s: out of memory allocating an array\n", "demo");
free((char *) its);
return (char **) 0;
}
for (i = 1; i < rows; ++i)
its[i] = &(its[0][i * cols * size]);
return its;
}
void save_halfcomplexmatrix(char *fname, int num_rows,
int num_cols, int col_pad, fftw_complex *mx)
{
int i, j;
FILE *f;
/* Open halfcomplexmatrix file.
*/
if (!(f = fopen(fname,"w"))) {
perror("Cannot open halfcomplexmatrix file.");
exit(1);
}
/* Save header
*/
if (fprintf(f,"%s\n", HALFCOMPLEXMATRIX_MAGIC) == EOF) {
perror("Cannot write halfcomplexmatrix file");
exit(1);
}
/* Save number of rows and cols.
*/
if (fprintf(f, "%d\n%d\n", num_rows, num_cols) == EOF) {
perror("Cannot write halfcomplexmatrix file");
exit(1);
}
/* write data.
*/
for (i = 0; i < num_rows; i++)
for (j = 0; j < num_cols; j++)
if (fprintf(f, "%e\t%e\n", mx[i*num_cols+j].re,
mx[i*num_cols+j].im) == EOF) {
perror("Cannot write halfcomplexmatrix file");
exit(1);
}
fclose(f);
}
fftw_real *extract_realsubimage(fftw_real *image, fftw_real *subimg,
int num_rows, int num_cols, int im_row_length,
int sub_row_length, int start_row,
int start_col)
{
fftw_real *subimage;
int row, col;
subimage = (subimg == NULL?
Alloc_dvector(num_rows * sub_row_length * sizeof(fftw_real)) :
subimg);
for (row = 0; row < num_rows; row++)
for (col = 0; col < num_cols; col++)
subimage[row*sub_row_length+col] =
image[(start_row+row)*im_row_length+start_col+col];
return(subimage);
}
void save_realmatrix(char *fname, int num_rows, int num_cols, int col_pad,
fftw_real *mx)
{
int i, j;
FILE *f;
/* Open realmatrix file.
*/
if (!(f = fopen(fname,"w"))) {
perror("Cannot open realmatrix file.");
exit(1);
}
/* Save header
*/
if (fprintf(f,"%s\n", REALMATRIX_MAGIC) == EOF) {
perror("Cannot write realmatrix file");
exit(1);
}
/* Save number of rows and cols.
*/
if (fprintf(f, "%d\n%d\n", num_rows, num_cols) == EOF) {
perror("Cannot write realmatrix file");
exit(1);
}
/* write data.
*/
for (i = 0; i < num_rows; i++)
for (j = 0; j < num_cols; j++)
if (fprintf(f, "%e\n", mx[i*(num_cols+col_pad)+j]) == EOF) {
perror("Cannot write realmatrix file");
exit(1);
}
fclose(f);
}
gray **real_to_gray(fftw_real *realim, gray **img, int num_rows, int num_cols,
int real_row_length)
{
double min, max, scale;
int row, col, k;
gray **im;
im = (img == NULL? pgm_allocarray(num_cols, num_rows) : img);
min = 100000000000;
max = -100000000000;
for (row = 0; row < num_rows; row++)
for (col = 0; col < num_cols; col++) {
k = row * real_row_length + col;
min = MIN(min, realim[k]);
max = MAX(max, realim[k]);
}
scale = 255.0 / (max - min);
for (row = 0; row < num_rows; row++)
for (col = 0; col < num_cols; col++)
im[row][col] = (int)(scale * (realim[row * real_row_length + col] - min));
return(im);
}
fftw_real *load_realmatrix(char *fname, int *num_rows, int *num_cols,
int col_pad)
{
int i, j;
FILE *f;
fftw_real *mx;
char buf[100];
/* Open realmatrix file.
*/
if (!(f = fopen(fname,"r"))) {
perror("Cannot open realmatrix file.");
exit(1);
}
/* Check realmatrix file type.
*/
if (fscanf(f,"%s\n", buf) != 1) {
perror("Cannot read from realmatrix file (wrong format?)");
exit(1);
}
if (strcmp(buf, REALMATRIX_MAGIC)) {
perror("Illegal realmatrix file!");
exit(1);
}
/* Read number of rows and cols.
*/
if (fscanf(f, "%d\n%d\n", num_rows, num_cols) == EOF) {
perror("Cannot read realmatrix file");
exit(1);
}
/* Allocate matrix and read data.
*/
mx = Alloc_dvector(*num_rows * (*num_cols+col_pad));
for (i = 0; i < *num_rows; i++)
for (j = 0; j < *num_cols; j++)
/* modify 1 */
if (fscanf(f, "%lf\n", mx+(i*(*num_cols+col_pad)+j)) == EOF) {
perror("Cannot read realmatrix file");
exit(1);
}
fclose(f);
return(mx);
}
void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with dvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void pm_freearray(char **its, int rows)
{
free(its[0]);
free(its);
}
double *dvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
double *v;
v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
// if (!v) nrerror("valami");
return v-nl+NR_END;
}
int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
int *v;
v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
// if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}
fftw_real *real_scaling(fftw_real *realim, fftw_real *out, int num_rows,
int num_cols, int real_row_length, int out_row_length)
{
double min, max, mean, scale;
int row, col, k;
fftw_real *scaled;
scaled = (out == NULL? Alloc_dvector(num_rows * out_row_length) : out);
/* Scale input image to [-1,1] and set the mean value to zero.
*/
min = 100000;
max = -100000;
mean = 0.0;
for (row = 0; row < num_rows; row++)
for (col = 0; col < num_cols; col++) {
k = row*real_row_length + col;
min = MIN(min, realim[k]);
max = MAX(max, realim[k]);
mean += realim[k];
}
mean /= ((double)num_cols * num_rows);
scale = MAX(max-mean, mean-min);
for (row = 0; row < num_rows; row++)
for (col = 0; col < num_cols; col++)
scaled[row*out_row_length + col] =
(realim[row*real_row_length + col] - mean) / scale;
return(scaled);
}
/**************************************
* mirror_pad_real()
**************************************/
fftw_real *mirror_pad_real(fftw_real *image, fftw_real *pad, int num_rows,
int num_cols, int im_row_length, int pad_row_length,
int row_pad, int col_pad)
{
fftw_real *padded;
int row, col;
padded = (pad == NULL? Alloc_dvector((num_rows+2*row_pad) *
pad_row_length * sizeof(fftw_real)) :
pad);
/* Copy image.
*/
for (row = 0; row < num_rows; row++)
for (col = 0; col < num_cols; col++)
padded[(row+row_pad)*pad_row_length + col+col_pad] =
image[row*im_row_length + col];
/* Mirror rows.
*/
for (row = 1; row <= row_pad; row++)
for (col = col_pad; col < num_cols + col_pad; col++) {
padded[(row_pad - row)*pad_row_length + col] =
padded[(row_pad + row)*pad_row_length + col];
padded[(row_pad + num_rows-1 + row)*pad_row_length + col] =
padded[(row_pad + num_rows-1 - row)*pad_row_length + col];
}
/* Mirror cols.
*/
for (row = 0; row < num_rows + 2*row_pad; row++)
for (col = 1; col <= col_pad; col++) {
padded[row*pad_row_length + col_pad-col] =
padded[row*pad_row_length + col_pad+col];
padded[row*pad_row_length + col_pad + num_cols-1 + col] =
padded[row*pad_row_length + col_pad + num_cols-1 - col];
}
return(padded);
}
void ft_even_gabor(fftw_complex *filter, int num_rows, int num_cols,
double U, double SigmaU, double SigmaV, double theta)
{
int i, j; /* Coordinates in the output frequency domain matrix*/
double u, v; /* Corresponding coords. in non-rotated image */
double sinTheta, cosTheta;
double halfInvSigmaUsq, halfInvSigmaVsq, A;
double iCont, jCont; /* Frequencies */
double t;
int ft_num_cols = num_cols/2+1; /*Num. columns of the half-complex output*/
/* Herein, we use the fact that the Fourier transform of a real and
* even function is also real and even. Thus, it is enogh to compute
* the coefficients in the positiv quarter of the frequency domain:
* [0,1/num_cols], [0,1/num_rows]. However, we will store the
* coefficients of the positiv half [0,1/num_cols],
* [-1/num_rows,1/num_rows] because FFTW uses the same format. Thus
* convolution and inverse FFT are easier to implement.
*/
//sincos(theta, &sinTheta, &cosTheta);
sinTheta = sin( theta );
cosTheta = cos( theta );
A = 1 / (2.0 * M_PI * SigmaU * SigmaV);
halfInvSigmaUsq = 0.5 / SQ(SigmaU);
halfInvSigmaVsq = 0.5 / SQ(SigmaV);
for (i = 0; i < num_rows/2+1; i++) {
iCont = (double)i/num_rows;
for (j = 0; j < ft_num_cols; j++) {
jCont = (double)j/num_cols;
u = iCont * cosTheta + jCont * sinTheta;
v = jCont * cosTheta - iCont * sinTheta;
t = SQ(v * num_rows/num_cols) * halfInvSigmaVsq;
filter[i*ft_num_cols+j].re =
A * (exp(-SQ(u - U) * halfInvSigmaUsq - t)
+ exp(-SQ(u + U) * halfInvSigmaUsq - t));
filter[i*ft_num_cols+j].im = 0.0;
if (i) {
filter[(num_rows-i)*ft_num_cols+j].re = filter[i*ft_num_cols+j].re;
filter[(num_rows-i)*ft_num_cols+j].im = 0.0;
}
}
}
/* In order to ensure that the filters do not respond to regions
* with constant intensity, we set the MTF of the filter at (0,0)
* to zero. As a result, the filtered image has a zero mean.
*/
filter[0].re = filter[0].im = 0.0;
}
/**************************************
* halfcomplex_to_gray()
**************************************/
gray **halfcomplex_to_gray(fftw_complex *hcim, gray **img,
int num_rows, int num_cols, int hcim_row_length)
{
double min, max, scale;
int row, col, k;
gray **im;
im = (img == NULL? pgm_allocarray(num_cols, num_rows) : img);
/* Compute min and max values.
*/
min = 1000000; //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
max = -1000000;
for (row = 0; row < num_rows; row++)
for (col = 0; col < num_cols/2+1; col++) {
k = row * hcim_row_length + col;
min = MIN(min, hcim[k].re);
max = MAX(max, hcim[k].re);
}
/* Store frequency domain magnitudes in PGM (gray) format.
*/
scale = 255.0 / (max - min);
for (row = 0; row < num_rows/2+1; row++)
for (col = 0; col < num_cols/2+1; col++) {
if (row+num_rows/2 < num_rows)
im[row+num_rows/2][num_cols/2-col] =
(int)(scale * (hcim[row*hcim_row_length+col].re - min));
if (row+num_rows/2 < num_rows && col+num_cols/2 < num_cols)
im[row+num_rows/2][col+num_cols/2] =
im[row+num_rows/2][num_cols/2-col];
}
for(row--; row < num_rows; row++)
for (col = 0; col < num_cols/2+1; col++) {
im[row - num_rows/2][num_cols/2-col] =
(int)(scale * (hcim[row*hcim_row_length+col].re - min));
if (col+num_cols/2 < num_cols)
im[row - num_rows/2][col+num_cols/2] =
im[row - num_rows/2][num_cols/2-col];
}
return(im);
}
/*****************************
* inverse()
*****************************/
double inverse(double **a, double **res, int n)
{
double d,*col, **tmp;
int i,j,*indx;
if (n==1) {
res[0][0] = 1.0 / a[0][0];
return(a[0][0]);
}
col = Alloc_dvector(n);
indx = Alloc_ivector(n);
/* Keep input matrix unchanged
*/
tmp = Alloc_dmatrix(n,n);
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
tmp[i][j] = a[i][j];
/* Decompose the matrix just once. This returns d as +/- 1. */
// ludcmp(tmp,n,indx,&d);
/* Find inverse by columns. */
for(j=0;j= 0)
for (j=ii;j<=i-1;j++)
sum -= a[i][j]*b[j];
else if (sum)
ii=i;
/* A nonzero element was encountered, so from now on we will have
* to do the sums in the loop above.
*/
b[i]=sum;
}
for (i=n-1;i>=0;i--) {
/* Now we do the backsubstitution, equation (2.3.7). */
sum=b[i];
for (j=i+1;j big)
big=temp;
if (big == 0.0)
/* No nonzero largest element. */
// nrerror("Singular matrix in routine ludcmp");
vv[i]=1.0/big; /* Save the scaling. */
}
for (j=0;j= big) {
/* Is the gure of merit for the pivot better than the best so far? */
big=dum;
imax=i;
}
}
/* Do we need to interchange rows? */
if (j != imax) {
/* Yes, do so... */
for (k=0;k= 1; k--)
{
x[k]=0;
for ( j = k+1; j <= n; j++ )
{
x[k]=x[k]-a[k][j]*x[j];
}
x[k]=(y[k]+x[k])/a[k][k];
}
}
/* A m trix inverz‚nek kisz mit sa */
void inverz(double** a,double** e, int n)
{
int i,j;
double *x, *y;
x = new double[n+1];
y = new double[n+1];
for ( j = 1; j <= n; j++ )
{
for ( i = 1; i <= n; i++ )
{
a[i][n+1]=e[i][j];
}
yvektor(y,a,n);
xvektor(x,y,a,n);
for ( i = 1; i <= n; i++ )
{
e[i][j]=x[i];
}
}
// matrixkiir(e,n);
}
/* A determin ns kisz mit sa */
double det(double** a,int n, int d)
{
int k;
double c;
c=a[1][1];
for ( k = 2; k <= n; k++ )
{
c=c*a[k][k];
}
return (d * c);
}