1 /* forward FFT
2 *
3 * Author: Nicos Dessipris
4 * Written on: 12/04/1990
5 * Modified on : 09/05/1990 to cope with float input
6 * Modified on : 08/03/1991 history removed
7 * Modified on : 03/04/1991 to cope with any input
8 *
9 * 28/6/95 JC
10 * - rewritten to use im_clip2f() rather than own code
11 * - memory leaks fixed
12 * 10/9/98 JC
13 * - frees memory more quickly
14 * 2/4/02 JC
15 * - fftw code added
16 * 13/7/02 JC
17 * - output Type set to IM_TYPE_FOURIER to help nip
18 * 27/2/03 JC
19 * - exploits real_to_complex() path in libfftw for real input (thanks
20 * Matt) for a 2x speed-up
21 * 17/11/03 JC
22 * - fix a segv for wider than high images in the real_to_complex() path
23 * (thanks Andrey)
24 * - fixes to real_to_complex() path to give the correct result for
25 * non-square images, including odd widths and heights
26 * 3/11/04
27 * - added fftw3 support
28 * 7/2/10
29 * - cleanups
30 * - gtkdoc
31 * 25/3/10
32 * - have a "t" image linked to out to keep the image alive for longer
33 * 27/1/12
34 * - better setting of interpretation
35 * - remove own fft fallback code
36 * - remove fftw2 path
37 * - reduce memuse
38 * 3/1/14
39 * - redone as a class
40 */
41
42 /*
43
44 This file is part of VIPS.
45
46 VIPS is free software; you can redistribute it and/or modify
47 it under the terms of the GNU Lesser General Public License as published by
48 the Free Software Foundation; either version 2 of the License, or
49 (at your option) any later version.
50
51 This program is distributed in the hope that it will be useful,
52 but WITHOUT ANY WARRANTY; without even the implied warranty of
53 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
54 GNU Lesser General Public License for more details.
55
56 You should have received a copy of the GNU Lesser General Public License
57 along with this program; if not, write to the Free Software
58 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
59 02110-1301 USA
60
61 */
62
63 /*
64
65 These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
66
67 */
68
69 #ifdef HAVE_CONFIG_H
70 #include <config.h>
71 #endif /*HAVE_CONFIG_H*/
72 #include <vips/intl.h>
73
74 #include <stdio.h>
75 #include <stdlib.h>
76
77 #include <vips/vips.h>
78 #include <vips/internal.h>
79 #include "pfreqfilt.h"
80
81 #ifdef HAVE_FFTW
82
83 #include <fftw3.h>
84
85 typedef struct _VipsFwfft {
86 VipsFreqfilt parent_instance;
87
88 } VipsFwfft;
89
90 typedef VipsFreqfiltClass VipsFwfftClass;
91
92 G_DEFINE_TYPE( VipsFwfft, vips_fwfft, VIPS_TYPE_FREQFILT );
93
94 /* Real to complex forward transform.
95 */
96 static int
rfwfft1(VipsObject * object,VipsImage * in,VipsImage ** out)97 rfwfft1( VipsObject *object, VipsImage *in, VipsImage **out )
98 {
99 VipsFwfft *fwfft = (VipsFwfft *) object;
100 VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
101 VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( fwfft );
102 const guint64 size = VIPS_IMAGE_N_PELS( in );
103 const int half_width = in->Xsize / 2 + 1;
104
105 double *half_complex;
106 double *planner_scratch;
107
108 fftw_plan plan;
109 double *buf, *q, *p;
110 int x, y;
111
112 if( vips_check_mono( class->nickname, in ) ||
113 vips_check_uncoded( class->nickname, in ) )
114 return( -1 );
115
116 /* Convert input to a real double membuffer.
117 */
118 t[1] = vips_image_new_memory();
119 if( vips_cast_double( in, &t[0], NULL ) ||
120 vips_image_write( t[0], t[1] ) )
121 return( -1 );
122
123 /* Make the plan for the transform. Yes, they really do use nx for
124 * height and ny for width. Use a separate scratch buffer for the
125 * planner, we can't overwrite real->data
126 */
127 if( !(planner_scratch = VIPS_ARRAY( fwfft,
128 VIPS_IMAGE_N_PELS( in ), double )) )
129 return( -1 );
130 if( !(half_complex = VIPS_ARRAY( fwfft,
131 in->Ysize * half_width * 2, double )) )
132 return( -1 );
133 if( !(plan = fftw_plan_dft_r2c_2d( in->Ysize, in->Xsize,
134 planner_scratch, (fftw_complex *) half_complex,
135 0 )) ) {
136 vips_error( class->nickname,
137 "%s", _( "unable to create transform plan" ) );
138 return( -1 );
139 }
140
141 fftw_execute_dft_r2c( plan,
142 (double *) t[1]->data, (fftw_complex *) half_complex );
143
144 fftw_destroy_plan( plan );
145
146 /* Write to out as another memory buffer.
147 */
148 *out = vips_image_new_memory();
149 if( vips_image_pipelinev( *out, VIPS_DEMAND_STYLE_ANY, in, NULL ) )
150 return( -1 );
151 (*out)->BandFmt = VIPS_FORMAT_DPCOMPLEX;
152 (*out)->Type = VIPS_INTERPRETATION_FOURIER;
153 if( !(buf = VIPS_ARRAY( fwfft, VIPS_IMAGE_N_PELS( *out ), double )) )
154 return( -1 );
155
156 /* Copy and normalise. The right half is the up/down and
157 * left/right flip of the left, but conjugated. Do the first
158 * row separately, then mirror around the centre row.
159 */
160 p = half_complex;
161 q = buf;
162
163 for( x = 0; x < half_width; x++ ) {
164 q[0] = p[0] / size;
165 q[1] = p[1] / size;
166 p += 2;
167 q += 2;
168 }
169
170 p = half_complex + ((in->Xsize + 1) / 2 - 1) * 2;
171
172 for( x = half_width; x < (*out)->Xsize; x++ ) {
173 q[0] = p[0] / size;
174 q[1] = -1.0 * p[1] / size;
175 p -= 2;
176 q += 2;
177 }
178
179 if( vips_image_write_line( *out, 0, (VipsPel *) buf ) )
180 return( -1 );
181
182 for( y = 1; y < (*out)->Ysize; y++ ) {
183 p = half_complex + y * half_width * 2;
184 q = buf;
185
186 for( x = 0; x < half_width; x++ ) {
187 q[0] = p[0] / size;
188 q[1] = p[1] / size;
189 p += 2;
190 q += 2;
191 }
192
193 /* Good grief.
194 */
195 p = half_complex + 2 *
196 (((*out)->Ysize - y + 1) * half_width - 2 +
197 (in->Xsize & 1));
198
199 for( x = half_width; x < (*out)->Xsize; x++ ) {
200 q[0] = p[0] / size;
201 q[1] = -1.0 * p[1] / size;
202 p -= 2;
203 q += 2;
204 }
205
206 if( vips_image_write_line( *out, y, (VipsPel *) buf ) )
207 return( -1 );
208 }
209
210 return( 0 );
211 }
212
213 /* Complex to complex forward transform.
214 */
215 static int
cfwfft1(VipsObject * object,VipsImage * in,VipsImage ** out)216 cfwfft1( VipsObject *object, VipsImage *in, VipsImage **out )
217 {
218 VipsFwfft *fwfft = (VipsFwfft *) object;
219 VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
220 VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( fwfft );
221
222 fftw_plan plan;
223 double *planner_scratch;
224 double *buf, *q, *p;
225 int x, y;
226
227 if( vips_check_mono( class->nickname, in ) ||
228 vips_check_uncoded( class->nickname, in ) )
229 return( -1 );
230
231 /* Convert input to a complex double membuffer.
232 */
233 t[1] = vips_image_new_memory();
234 if( vips_cast_dpcomplex( in, &t[0], NULL ) ||
235 vips_image_write( t[0], t[1] ) )
236 return( -1 );
237
238 /* We have to have a separate buffer for the planner to work on.
239 */
240 if( !(planner_scratch = VIPS_ARRAY( fwfft,
241 VIPS_IMAGE_N_PELS( in ) * 2, double )) )
242 return( -1 );
243
244 /* Make the plan for the transform.
245 */
246 if( !(plan = fftw_plan_dft_2d( in->Ysize, in->Xsize,
247 (fftw_complex *) planner_scratch,
248 (fftw_complex *) planner_scratch,
249 FFTW_FORWARD,
250 0 )) ) {
251 vips_error( class->nickname,
252 "%s", _( "unable to create transform plan" ) );
253 return( -1 );
254 }
255
256 fftw_execute_dft( plan,
257 (fftw_complex *) t[1]->data, (fftw_complex *) t[1]->data );
258
259 fftw_destroy_plan( plan );
260
261 /* Write to out as another memory buffer.
262 */
263 *out = vips_image_new_memory();
264 if( vips_image_pipelinev( *out, VIPS_DEMAND_STYLE_ANY, in, NULL ) )
265 return( -1 );
266 (*out)->BandFmt = VIPS_FORMAT_DPCOMPLEX;
267 (*out)->Type = VIPS_INTERPRETATION_FOURIER;
268 if( !(buf = VIPS_ARRAY( fwfft, VIPS_IMAGE_N_PELS( *out ), double )) )
269 return( -1 );
270
271 /* Copy to out, normalise.
272 */
273 p = (double *) t[1]->data;
274 for( y = 0; y < (*out)->Ysize; y++ ) {
275 guint64 size = VIPS_IMAGE_N_PELS( *out );
276
277 q = buf;
278
279 for( x = 0; x < (*out)->Xsize; x++ ) {
280 q[0] = p[0] / size;
281 q[1] = p[1] / size;
282 p += 2;
283 q += 2;
284 }
285
286 if( vips_image_write_line( *out, y, (VipsPel *) buf ) )
287 return( -1 );
288 }
289
290 return( 0 );
291 }
292
293 static int
vips_fwfft_build(VipsObject * object)294 vips_fwfft_build( VipsObject *object )
295 {
296 VipsFreqfilt *freqfilt = VIPS_FREQFILT( object );
297 VipsFwfft *fwfft = (VipsFwfft *) object;
298 VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
299
300 VipsImage *in;
301
302 if( VIPS_OBJECT_CLASS( vips_fwfft_parent_class )->
303 build( object ) )
304 return( -1 );
305
306 in = freqfilt->in;
307
308 if( vips_image_decode( in, &t[0] ) )
309 return( -1 );
310 in = t[0];
311
312 if( vips_band_format_iscomplex( in->BandFmt ) ) {
313 if( vips__fftproc( VIPS_OBJECT( fwfft ), in, &t[1],
314 cfwfft1 ) )
315 return( -1 );
316 }
317 else {
318 if( vips__fftproc( VIPS_OBJECT( fwfft ), in, &t[1],
319 rfwfft1 ) )
320 return( -1 );
321 }
322
323 if( vips_image_write( t[1], freqfilt->out ) )
324 return( -1 );
325
326 return( 0 );
327 }
328
329 static void
vips_fwfft_class_init(VipsFwfftClass * class)330 vips_fwfft_class_init( VipsFwfftClass *class )
331 {
332 VipsObjectClass *vobject_class = VIPS_OBJECT_CLASS( class );
333
334 vobject_class->nickname = "fwfft";
335 vobject_class->description = _( "forward FFT" );
336 vobject_class->build = vips_fwfft_build;
337
338 }
339
340 static void
vips_fwfft_init(VipsFwfft * fwfft)341 vips_fwfft_init( VipsFwfft *fwfft )
342 {
343 }
344
345 #endif /*HAVE_FFTW*/
346
347 /**
348 * vips_fwfft: (method)
349 * @in: input image
350 * @out: (out): output image
351 * @...: %NULL-terminated list of optional named arguments
352 *
353 * Transform an image to Fourier space.
354 *
355 * VIPS uses the fftw Fourier Transform library. If this library was not
356 * available when VIPS was configured, these functions will fail.
357 *
358 * See also: vips_invfft().
359 *
360 * Returns: 0 on success, -1 on error.
361 */
362 int
vips_fwfft(VipsImage * in,VipsImage ** out,...)363 vips_fwfft( VipsImage *in, VipsImage **out, ... )
364 {
365 va_list ap;
366 int result;
367
368 va_start( ap, out );
369 result = vips_call_split( "fwfft", ap, in, out );
370 va_end( ap );
371
372 return( result );
373 }
374
375