1 /**********
2 *
3 * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
4 * University of Manchester, UK. All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without modification,
7 * are permitted provided that the following conditions are met:
8 *
9 * . Redistributions of source code must retain the above copyright notice,
10 * this list of conditions and the following disclaimer.
11 *
12 * . Redistributions in binary form must reproduce the above copyright notice,
13 * this list of conditions and the following disclaimer in the documentation
14 * and/or other materials provided with the distribution.
15 *
16 * . Neither the name of the University of Manchester nor the names of its
17 * contributors may be used to endorse or promote products derived from this
18 * software without specific prior written permission.
19 *
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 *
33 **********
34 *
35 * Program : TINA
36 * File : $Source: /home/tina/cvs/tina-libs/tina/image/imgPrc_gabor.c,v $
37 * Date : $Date: 2003/09/30 16:53:49 $
38 * Version : $Revision: 1.6 $
39 * CVS Id : $Id: imgPrc_gabor.c,v 1.6 2003/09/30 16:53:49 ipoole Exp $
40 *
41 * Author : Legacy TINA
42 */
43
44 /**
45 * @file
46 * @brief Gabor filter and image convolution.
47 *
48 * Gabor filters are defined in the Fourier domain by harmonic functions modulated
49 * by a Gaussian distribution.
50 *
51 */
52
53 #include "imgPrc_gabor.h"
54
55 #if HAVE_CONFIG_H
56 #include <config.h>
57 #endif
58
59 #include <math.h>
60 #include <tina/sys/sysDef.h>
61 #include <tina/sys/sysPro.h>
62 #include <tina/math/mathDef.h>
63 #include <tina/math/mathPro.h>
64 #include <tina/image/img_GenDef.h>
65 #include <tina/image/img_GenPro.h>
66 #include <tina/image/imgPrc_conv_1d.h>
67 #include <tina/image/imgPrc_prof1.h>
68 #include <tina/image/imgPrc_fourier.h>
69
70 Prof1 *gabor_profile(float phase, float sigma, float omega, float nsigma)
71 {
72 float *el;
73 Prof1 *prof;
74 float k;
75 int i, n;
76
77 if (nsigma <= 0.0)
78 return (NULL);
79
80 /** allocate nsigma*sigma mask **/
81 n = (int) (nsigma * sigma);
82 prof = prof1_alloc(-n, n + 1, float_v);
83 el = prof->el.float_v;
84
85 if (n == 0)
86 {
87 /** return delta function **/
88 el[0] = (float) 1.0;
89 return (prof);
90 }
91 k = (float) (1.0 / (2.0 * sigma * sigma));
92 for (i = -n; i < n + 1; i++)
93 el[i] = (float) (exp(-k * i * i) * sin(omega * i + phase));
94
95 return (prof);
96 }
97
98 Imrect *im_gabor(Imrect * im1, double phasex, double sigmax, double omegax, double nsigmax, double phasey, double sigmay, double omegay, double nsigmay)
99 {
100 Prof1 *profx = gabor_profile((float) phasex, (float) sigmax, (float) omegax, (float) nsigmax);
101 Prof1 *profy = gabor_profile((float) phasey, (float) sigmay, (float) omegay, (float) nsigmay);
102 Imrect *im2;
103
104 im2 = im_conv_separable(im1, profx, profy);
105 prof1_free(profx);
106 prof1_free(profy);
107 return (im2);
108 }
109
110 /*
111 * applies a convolution of 2D gabor filter: k = centre frequency b = 1-sigma bandwidth
112 * in octaves theta = orientation, to im.
113 */
114
115 Imrect *im_gabor_fft(Imrect * im, double k, double b, double theta)
116 {
117 Imrect *im1;
118 Imrect *im2;
119 double s, kx, ky, c1, c2;
120 double *px, *py;
121 int x, y, wx, wy;
122
123 /* work in frequency space */
124 im1 = im_fft(im, im->region);
125 wx = im->width;
126 wy = im->height;
127 b = pow(2.0, b);
128 s = k * (b - 1.0) / (b + 1.0);
129 kx = k * cos(theta);
130 ky = k * sin(theta);
131 c1 = 1.0 / (s * sqrt(TWOPI));
132 c2 = -0.5 / (s * s);
133
134 /** set up frequency profiles in x & y (with wrap-around) */
135 px = tvector_alloc(0, wx, double);
136 for (x = 0; x < wx; x++)
137 {
138 int x1, x2;
139
140 x1 = (int) fabs((double) x - kx);
141 x2 = wx - x1;
142 x1 = MIN(x1, x2);
143 px[x] = c1 * exp(c2 * x1 * x1);
144 }
145 py = tvector_alloc(0, wy, double);
146 for (y = 0; y < wy; y++)
147 {
148 int y1, y2;
149
150 y1 = (int) fabs((double) y - ky);
151 y2 = wy - y1;
152 y1 = MIN(y1, y2);
153 py[y] = c1 * exp(c2 * y1 * y1);
154 }
155
156 /* apply convolutuion as product in frequency space */
157 for (x = 0; x < wx; x++)
158 for (y = 0; y < wy; y++)
159 {
160 Complex z = {Complex_id};
161
162 IM_PIX_GETZ(im1, y, x, z);
163 z = cmplx_times(px[x] * py[y], z);
164 IM_PIX_SETZ(im1, y, x, z);
165 }
166 tvector_free(px, 0, double);
167 tvector_free(py, 0, double);
168 im2 = im_fft_inverse(im1, im1->region);
169 im_free(im1);
170 return (im2);
171 }
172
173 Imrect *im_fgabor(Imregion * roi, double k, double b, double theta)
174 {
175 Imrect *im1;
176 double s, kx, ky, c1, c2;
177 double *px, *py;
178 int x, y, wx, wy;
179 Complex z = {Complex_id};
180
181 /* work in frequency space */
182 wx = roi->ux - roi->lx;
183 wy = roi->uy - roi->ly;
184 im1 = im_alloc(wy, wx, NULL, complex_v);
185
186 b = pow(2.0, b);
187 s = k * (b - 1.0) / (b + 1.0);
188 kx = k * cos(theta);
189 ky = k * sin(theta);
190 c1 = 1.0 / (s * sqrt(TWOPI));
191 c2 = -0.5 / (s * s);
192
193 /** set up frequency profiles in x & y (with wrap-around) */
194 px = tvector_alloc(0, wx, double);
195 for (x = 0; x < wx; x++)
196 {
197 int x1, x2;
198
199 x1 = (int) fabs((double) x - kx);
200 x2 = wx - x1;
201 x1 = MIN(x1, x2);
202 px[x] = c1 * exp(c2 * x1 * x1);
203 }
204 py = tvector_alloc(0, wy, double);
205 for (y = 0; y < wy; y++)
206 {
207 int y1, y2;
208
209 y1 = (int) fabs((double) y - ky);
210 y2 = wy - y1;
211 y1 = MIN(y1, y2);
212 py[y] = c1 * exp(c2 * y1 * y1);
213 }
214
215 /* apply convolutuion as product in frequency space */
216 for (x = 0; x < wx; x++)
217 for (y = 0; y < wy; y++)
218 {
219 z = cmplx_unit();
220 z = cmplx_times(px[x] * py[y], z);
221 IM_PIX_SETZ(im1, y, x, z);
222 }
223 tvector_free(px, 0, double);
224 tvector_free(py, 0, double);
225 return (im1);
226 }
227
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.