/* Lear's GIST implementation, version 1.1, (c) INRIA 2009, Licence: PSFL */
/*--------------------------------------------------------------------------*/
#ifdef USE_GIST
/*--------------------------------------------------------------------------*/
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <fftw3.h>
#include <pthread.h>
#ifdef STANDALONE_GIST
#include "gist.h"
#else
#include <image.h>
#include <descriptors.h>
#endif
/*--------------------------------------------------------------------------*/
static pthread_mutex_t fftw_mutex = PTHREAD_MUTEX_INITIALIZER;
static void fftw_lock(void)
{
pthread_mutex_lock(&fftw_mutex);
}
static void fftw_unlock(void)
{
pthread_mutex_unlock(&fftw_mutex);
}
/*--------------------------------------------------------------------------*/
static image_t *image_add_padding(image_t *src, int padding)
{
int i, j;
image_t *img = image_new(src->width + 2*padding, src->height + 2*padding);
memset(img->data, 0x00, img->stride*img->height*sizeof(float));
for(j = 0; j < src->height; j++)
{
for(i = 0; i < src->width; i++) {
img->data[(j+padding)*img->stride+i+padding] = src->data[j*src->stride+i];
}
}
for(j = 0; j < padding; j++)
{
for(i = 0; i < src->width; i++)
{
img->data[j*img->stride+i+padding] = src->data[(padding-j-1)*src->stride+i];
img->data[(j+padding+src->height)*img->stride+i+padding] = src->data[(src->height-j-1)*src->stride+i];
}
}
for(j = 0; j < img->height; j++)
{
for(i = 0; i < padding; i++)
{
img->data[j*img->stride+i] = img->data[j*img->stride+padding+padding-i-1];
img->data[j*img->stride+i+padding+src->width] = img->data[j*img->stride+img->width-padding-i-1];
}
}
return img;
}
/*--------------------------------------------------------------------------*/
static color_image_t *color_image_add_padding(color_image_t *src, int padding)
{
int i, j;
color_image_t *img = color_image_new(src->width + 2*padding, src->height + 2*padding);
for(j = 0; j < src->height; j++)
{
for(i = 0; i < src->width; i++)
{
img->c1[(j+padding)*img->width+i+padding] = src->c1[j*src->width+i];
img->c2[(j+padding)*img->width+i+padding] = src->c2[j*src->width+i];
img->c3[(j+padding)*img->width+i+padding] = src->c3[j*src->width+i];
}
}
for(j = 0; j < padding; j++)
{
for(i = 0; i < src->width; i++)
{
img->c1[j*img->width+i+padding] = src->c1[(padding-j-1)*src->width+i];
img->c2[j*img->width+i+padding] = src->c2[(padding-j-1)*src->width+i];
img->c3[j*img->width+i+padding] = src->c3[(padding-j-1)*src->width+i];
img->c1[(j+padding+src->height)*img->width+i+padding] = src->c1[(src->height-j-1)*src->width+i];
img->c2[(j+padding+src->height)*img->width+i+padding] = src->c2[(src->height-j-1)*src->width+i];
img->c3[(j+padding+src->height)*img->width+i+padding] = src->c3[(src->height-j-1)*src->width+i];
}
}
for(j = 0; j < img->height; j++)
{
for(i = 0; i < padding; i++)
{
img->c1[j*img->width+i] = img->c1[j*img->width+padding+padding-i-1];
img->c2[j*img->width+i] = img->c2[j*img->width+padding+padding-i-1];
img->c3[j*img->width+i] = img->c3[j*img->width+padding+padding-i-1];
img->c1[j*img->width+i+padding+src->width] = img->c1[j*img->width+img->width-padding-i-1];
img->c2[j*img->width+i+padding+src->width] = img->c2[j*img->width+img->width-padding-i-1];
img->c3[j*img->width+i+padding+src->width] = img->c3[j*img->width+img->width-padding-i-1];
}
}
return img;
}
/*--------------------------------------------------------------------------*/
static void image_rem_padding(image_t *dest, image_t *src, int padding)
{
int i, j;
for(j = 0; j < dest->height; j++)
{
for(i = 0; i < dest->width; i++) {
dest->data[j*dest->stride+i] = src->data[(j+padding)*src->stride+i+padding];
}
}
}
/*--------------------------------------------------------------------------*/
static void color_image_rem_padding(color_image_t *dest, color_image_t *src, int padding)
{
int i, j;
for(j = 0; j < dest->height; j++)
{
for(i = 0; i < dest->width; i++)
{
dest->c1[j*dest->width+i] = src->c1[(j+padding)*src->width+i+padding];
dest->c2[j*dest->width+i] = src->c2[(j+padding)*src->width+i+padding];
dest->c3[j*dest->width+i] = src->c3[(j+padding)*src->width+i+padding];
}
}
}
/*--------------------------------------------------------------------------*/
static void fftshift(float *data, int w, int h)
{
int i, j;
float *buff = (float *) malloc(w*h*sizeof(float));
memcpy(buff, data, w*h*sizeof(float));
for(j = 0; j < (h+1)/2; j++)
{
for(i = 0; i < (w+1)/2; i++) {
data[(j+h/2)*w + i+w/2] = buff[j*w + i];
}
for(i = 0; i < w/2; i++) {
data[(j+h/2)*w + i] = buff[j*w + i+(w+1)/2];
}
}
for(j = 0; j < h/2; j++)
{
for(i = 0; i < (w+1)/2; i++) {
data[j*w + i+w/2] = buff[(j+(h+1)/2)*w + i];
}
for(i = 0; i < w/2; i++) {
data[j*w + i] = buff[(j+(h+1)/2)*w + i+(w+1)/2];
}
}
free(buff);
}
/*--------------------------------------------------------------------------*/
static image_list_t *create_gabor(int nscales, const int *or, int width, int height)
{
int i, j, fn;
image_list_t *G = image_list_new();
int nfilters = 0;
for(i=0;i<nscales;i++) nfilters+=or[i];
float **param = (float **) malloc(nscales * nfilters * sizeof(float *));
for(i = 0; i < nscales * nfilters; i++) {
param[i] = (float *) malloc(4*sizeof(float));
}
float *fx = (float *) malloc(width*height*sizeof(float));
float *fy = (float *) malloc(width*height*sizeof(float));
float *fr = (float *) malloc(width*height*sizeof(float));
float *f = (float *) malloc(width*height*sizeof(float));
int l = 0;
for(i = 1; i <= nscales; i++)
{
for(j = 1; j <= or[i-1]; j++)
{
param[l][0] = 0.35f;
param[l][1] = 0.3/pow(1.85f, i-1);
param[l][2] = 16*pow(or[i-1], 2)/pow(32, 2);
param[l][3] = M_PI/(or[i-1])*(j-1);
l++;
}
}
for(j = 0; j < height; j++)
{
for(i = 0; i < width; i++)
{
fx[j*width + i] = (float) i - width/2.0f;
fy[j*width + i] = (float) j - height/2.0f;
fr[j*width + i] = sqrt(fx[j*width + i]*fx[j*width + i] + fy[j*width + i]*fy[j*width + i]);
f[j*width + i] = atan2(fy[j*width + i], fx[j*width + i]);
}
}
fftshift(fr, width, height);
fftshift(f, width, height);
for(fn = 0; fn < nfilters; fn++)
{
image_t *G0 = image_new(width, height);
float *f_ptr = f;
float *fr_ptr = fr;
for(j = 0; j < height; j++)
{
for(i = 0; i < width; i++)
{
float tmp = *f_ptr++ + param[fn][3];
if(tmp < -M_PI) {
tmp += 2.0f*M_PI;
}
else if (tmp > M_PI) {
tmp -= 2.0f*M_PI;
}
G0->data[j*G0->stride+i] = exp(-10.0f*param[fn][0]*(*fr_ptr/height/param[fn][1]-1)*(*fr_ptr/width/param[fn][1]-1)-2.0f*param[fn][2]*M_PI*tmp*tmp);
fr_ptr++;
}
}
image_list_append(G, G0);
}
for(i = 0; i < nscales * nfilters; i++) {
free(param[i]);
}
free(param);
free(fx);
free(fy);
free(fr);
free(f);
return G;
}
/*------------