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
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.