~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/vision/improc/im_gabor.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**@(#)
  2 **/
  3 /* gabor.c
  4  * 
  5  * gabor filter functions
  6  * 
  7  */
  8 
  9 #include <math.h>
 10 #include <tina/sys.h>
 11 #include <tina/sysfuncs.h>
 12 #include <tina/math.h>
 13 #include <tina/mathfuncs.h>
 14 #include <tina/vision.h>
 15 #include <tina/visionfuncs.h>
 16 
 17 Prof1  *prof1_alloc();
 18 Imrect *im_conv_separable();
 19 
 20 Prof1  *gabor_profile(float phase, float sigma, float omega, float nsigma)
 21 {
 22     float  *el;
 23     Prof1  *prof;
 24     float   k;
 25     int     i, n;
 26 
 27     if (nsigma <= 0.0)
 28         return (NULL);
 29 
 30     /** allocate nsigma*sigma mask **/
 31     n = (int) (nsigma * sigma);
 32     prof = prof1_alloc(-n, n + 1, float_v);
 33     el = prof->el.float_v;
 34 
 35     if (n == 0)
 36     {
 37         /** return delta function **/
 38         el[0] = (float) 1.0;
 39         return (prof);
 40     }
 41     k = (float) (1.0 / (2.0 * sigma * sigma));
 42     for (i = -n; i < n + 1; i++)
 43         el[i] = (float) (exp(-k * i * i) * sin(omega * i + phase));
 44 
 45     return (prof);
 46 }
 47 
 48 Imrect *im_gabor(Imrect * im1, double phasex, double sigmax, double omegax, double nsigmax, double phasey, double sigmay, double omegay, double nsigmay)
 49 {
 50     Prof1  *profx = gabor_profile((float)phasex, (float)sigmax, (float)omegax, (float)nsigmax);
 51     Prof1  *profy = gabor_profile((float)phasey, (float)sigmay, (float)omegay, (float)nsigmay);
 52     Imrect *im2;
 53 
 54     im2 = im_conv_separable(im1, profx, profy);
 55     prof1_free(profx);
 56     prof1_free(profy);
 57     return (im2);
 58 }
 59 
 60 /* fft version of 2D gabor filter: k = centre frequency b = 1-sigma
 61  * bandwidth in octaves theta = orientation */
 62 
 63 Imrect *im_gabor_fft(Imrect * im, double k, double b, double theta)
 64 {
 65     Imrect *im1;
 66     Imrect *im2;
 67     double  s, kx, ky, c1, c2;
 68     double *px, *py;
 69     int     x, y, wx, wy;
 70 
 71     /* work in frequency space */
 72     im1 = im_fft(im,im->region);
 73     wx = im->width;
 74     wy = im->height;
 75     b = pow(2.0, b);
 76     s = k * (b - 1.0) / (b + 1.0);
 77     kx = k * cos(theta);
 78     ky = k * sin(theta);
 79     c1 = 1.0 / (s * sqrt(TWOPI));
 80     c2 = -0.5 / (s * s);
 81 
 82     /** set up frequency profiles in x & y (with wrap-around) */
 83     px = tvector_alloc(0, wx, double);
 84     for (x = 0; x < wx; x++)
 85     {
 86         int     x1, x2;
 87 
 88         x1 = (int)fabs((double)x - kx);
 89         x2 = wx - x1;
 90         x1 = MIN(x1, x2);
 91         px[x] = c1 * exp(c2 * x1 * x1);
 92     }
 93     py = tvector_alloc(0, wy, double);
 94     for (y = 0; y < wy; y++)
 95     {
 96         int     y1, y2;
 97 
 98         y1 = (int)fabs((double)y - ky);
 99         y2 = wy - y1;
100         y1 = MIN(y1, y2);
101         py[y] = c1 * exp(c2 * y1 * y1);
102     }
103 
104     /* apply convolutuion as product in frequency space */
105     for (x = 0; x < wx; x++)
106         for (y = 0; y < wy; y++)
107         {
108             Complex z = {Complex_id};
109 
110             IM_PIX_GETZ(im1, y, x, z);
111             z = cmplx_times(px[x] * py[y], z);
112             IM_PIX_SETZ(im1, y, x, z);
113         }
114     tvector_free(px, 0, double);
115     tvector_free(py, 0, double);
116     im2 = im_fft_inverse(im1,im1->region);
117     im_free(im1);
118     return (im2);
119 }
120 
121 Imrect *im_fgabor(Imregion *roi, double k, double b, double theta)
122 {
123     Imrect *im1;
124     double  s, kx, ky, c1, c2;
125     double *px, *py;
126     int     x, y, wx, wy;
127     Complex z = {Complex_id};
128    
129     /* work in frequency space */
130     wx = roi->ux - roi->lx;
131     wy = roi->uy - roi->ly;
132     im1 = im_alloc(wy,wx,NULL,complex_v);
133 
134     b = pow(2.0, b);
135     s = k * (b - 1.0) / (b + 1.0);
136     kx = k * cos(theta);
137     ky = k * sin(theta);
138     c1 = 1.0 / (s * sqrt(TWOPI));
139     c2 = -0.5 / (s * s);
140 
141     /** set up frequency profiles in x & y (with wrap-around) */
142     px = tvector_alloc(0, wx, double);
143     for (x = 0; x < wx; x++)
144     {
145         int     x1, x2;
146 
147         x1 = (int)fabs((double)x - kx);
148         x2 = wx - x1;
149         x1 = MIN(x1, x2);
150         px[x] = c1 * exp(c2 * x1 * x1);
151     }
152     py = tvector_alloc(0, wy, double);
153     for (y = 0; y < wy; y++)
154     {
155         int     y1, y2;
156 
157         y1 = (int)fabs((double)y - ky);
158         y2 = wy - y1;
159         y1 = MIN(y1, y2);
160         py[y] = c1 * exp(c2 * y1 * y1);
161     }
162 
163     /* apply convolutuion as product in frequency space */
164     for (x = 0; x < wx; x++)
165         for (y = 0; y < wy; y++)
166         {
167             z = cmplx_unit();
168             z = cmplx_times(px[x] * py[y], z);
169             IM_PIX_SETZ(im1, y, x, z);
170         }
171     tvector_free(px, 0, double);
172     tvector_free(py, 0, double);
173     return (im1);
174 }
175 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.