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

Linux Cross Reference
Tina6/tina-libs/tina/image/imgPrc_gabor.c

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

  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 

~ [ 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.