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

Linux Cross Reference
Tina6/tina-libs/tina/image/imgPrc_fourier.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_fourier.c,v $
 37  * Date    :  $Date: 2003/09/22 16:09:02 $
 38  * Version :  $Revision: 1.4 $
 39  * CVS Id  :  $Id: imgPrc_fourier.c,v 1.4 2003/09/22 16:09:02 tony Exp $
 40  *
 41  * Author  : Legacy TINA
 42  */
 43 
 44 /** 
 45  *  @file
 46  *  @brief Fast Fourier Transform (FFT) and inverse.
 47  *
 48  */
 49 
 50 #include "imgPrc_fourier.h"
 51 
 52 #if HAVE_CONFIG_H
 53 #include <config.h>
 54 #endif
 55 
 56 #include <math.h>
 57 #include <tina/sys/sysDef.h>
 58 #include <tina/sys/sysPro.h>
 59 #include <tina/math/mathDef.h>
 60 #include <tina/math/mathPro.h>
 61 #include <tina/image/imgDef.h>
 62 #include <tina/image/imgPro.h>
 63 
 64 Imrect         *im_fft(Imrect * im1, Imregion * region)
 65 {
 66         Imrect         *im2;
 67         Imregion        roi;
 68         Complex        *line;
 69         int             nx, ny;
 70         int             lx, ux, ly, uy;
 71         int             i;
 72 
 73         if (im1 == NULL)
 74                 return (NULL);
 75 
 76         if (region == NULL)
 77                 region = im1->region;
 78         if (region == NULL)
 79                 return (NULL);
 80 
 81         lx = roi.lx = region->lx;
 82         ux = roi.ux = region->ux;
 83         ly = roi.ly = region->ly;
 84         uy = roi.uy = region->uy;
 85 
 86         for (nx = 1; nx <= ux - lx; nx *= 2);
 87         for (ny = 1; ny <= uy - ly; ny *= 2);
 88         nx /= 2;
 89         ny /= 2;
 90 
 91         ux = roi.ux = lx + nx;
 92         uy = roi.uy = ly + ny;
 93 
 94         im2 = im_alloc(ny, nx, &roi, complex_v);
 95         line = zvector_alloc(lx, ux);
 96 
 97         for (i = ly; i < uy; ++i)
 98         {
 99                 im_get_rowz(line, im1, i, lx, ux);
100                 line += lx;
101                 fft_cmplx_inplace(line, nx);
102                 line -= lx;
103                 im_put_rowz(line, im2, i, lx, ux);
104         }
105 
106         zvector_free(line, lx);
107         line = zvector_alloc(ly, uy);
108         for (i = lx; i < ux; ++i)
109         {
110                 im_get_colz(line, im2, i, ly, uy);
111                 line += ly;
112                 fft_cmplx_inplace(line, ny);
113                 line -= ly;
114                 im_put_colz(line, im2, i, ly, uy);
115         }
116 
117         zvector_free(line, ly);
118         return (im2);
119 }
120 
121 Imrect         *im_fft_inverse(Imrect * im1, Imregion * region)
122 {
123         Imrect         *im2;
124         Imregion        roi;
125         Complex        *line;
126         int             nx, ny;
127         int             lx, ux, ly, uy;
128         int             i;
129 
130         if (im1 == NULL)
131                 return (NULL);
132 
133         if (region == NULL)
134                 region = im1->region;
135         if (region == NULL)
136                 return (NULL);
137 
138         lx = roi.lx = region->lx;
139         ux = roi.ux = region->ux;
140         ly = roi.ly = region->ly;
141         uy = roi.uy = region->uy;
142 
143         for (nx = 1; nx <= ux - lx; nx *= 2);
144         for (ny = 1; ny <= uy - ly; ny *= 2);
145         nx /= 2;
146         ny /= 2;
147 
148         ux = roi.ux = lx + nx;
149         uy = roi.uy = ly + ny;
150 
151         im2 = im_alloc(ny, nx, &roi, complex_v);
152         line = zvector_alloc(lx, ux);
153 
154         for (i = ly; i < uy; ++i)
155         {
156                 im_get_rowz(line, im1, i, lx, ux);
157                 line += lx;
158                 fft_inverse_cmplx_inplace(line, nx);
159                 line -= lx;
160                 im_put_rowz(line, im2, i, lx, ux);
161         }
162 
163         zvector_free(line, lx);
164         line = zvector_alloc(ly, uy);
165         for (i = lx; i < ux; ++i)
166         {
167                 im_get_colz(line, im2, i, ly, uy);
168                 line += ly;
169                 fft_inverse_cmplx_inplace(line, ny);
170                 line -= ly;
171                 im_put_colz(line, im2, i, ly, uy);
172         }
173 
174         zvector_free(line, ly);
175         return (im2);
176 }
177 
178 static void     line_power_spectrum(Complex * l1, float *l2, int from, int to)
179 {
180         int             i;
181 
182         for (i = from; i < to; i++)
183                 l2[i] = (float) cmplx_sqrmod(l1[i]);
184 }
185 
186 Imrect         *im_power_spectrum(Imrect * im1)
187 {
188         Imrect         *im2;
189         Complex        *l1;
190         float          *l2;
191         int             i;
192         int             ux, uy, hx, hy;
193 
194         if (im1 == NULL)
195                 return (NULL);
196 
197         ux = im1->width;
198         hx = ux / 2;
199         uy = im1->height;
200         hy = uy / 2;
201 
202         im2 = im_alloc(uy, ux, NULL, float_v);
203         l1 = zvector_alloc(0, ux);
204         l2 = fvector_alloc(0, ux);
205 
206         for (i = 0; i < hy; ++i)
207         {
208                 im_get_rowz(l1, im1, i, 0, hx);
209                 line_power_spectrum(l1, l2, 0, hx);
210                 im_put_rowf(l2 - hx, im2, i + hy, hx, ux);
211 
212                 im_get_rowz(l1, im1, i, hx, ux);
213                 line_power_spectrum(l1, l2, hx, ux);
214                 im_put_rowf(l2 + hx, im2, i + hy, 0, hx);
215         }
216 
217         for (i = hy; i < uy; ++i)
218         {
219                 im_get_rowz(l1, im1, i, 0, hx);
220                 line_power_spectrum(l1, l2, 0, hx);
221                 im_put_rowf(l2 - hx, im2, i - hy, hx, ux);
222 
223                 im_get_rowz(l1, im1, i, hx, ux);
224                 line_power_spectrum(l1, l2, hx, ux);
225                 im_put_rowf(l2 + hx, im2, i - hy, 0, hx);
226         }
227 
228         zvector_free(l1, 0);
229         fvector_free(l2, 0);
230 
231         return (im2);
232 }
233 

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