1 /* Inverse FFT
2 *
3 * Author: Nicos Dessipris
4 * Written on: 12/04/1990
5 * Modified on :
6 * 28/6/95 JC
7 * - rewritten, based on new im_invfft() code
8 * 10/9/98 JC
9 * - frees memory more quickly
10 * 2/4/02 JC
11 * - fftw code added
12 * 13/7/02 JC
13 * - Type reset
14 * 27/2/03 JC
15 * - tiny speed-up ... save 1 copy on write
16 * 22/1/04 JC
17 * - oops, fix for segv on wider than high fftw transforms
18 * 3/11/04
19 * - added fftw3 support
20 * 7/2/10
21 * - gtkdoc
22 * 27/1/12
23 * - better setting of interpretation
24 * - remove own fft fallback code
25 * - remove fftw2 path
26 * - reduce memuse
27 * 3/1/14
28 * - redone as a class
29 */
30
31 /*
32
33 This file is part of VIPS.
34
35 VIPS is free software; you can redistribute it and/or modify
36 it under the terms of the GNU Lesser General Public License as published by
37 the Free Software Foundation; either version 2 of the License, or
38 (at your option) any later version.
39
40 This program is distributed in the hope that it will be useful,
41 but WITHOUT ANY WARRANTY; without even the implied warranty of
42 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
43 GNU Lesser General Public License for more details.
44
45 You should have received a copy of the GNU Lesser General Public License
46 along with this program; if not, write to the Free Software
47 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
48 02110-1301 USA
49
50 */
51
52 /*
53
54 These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
55
56 */
57
58 #ifdef HAVE_CONFIG_H
59 #include <config.h>
60 #endif /*HAVE_CONFIG_H*/
61 #include <vips/intl.h>
62
63 #include <stdio.h>
64 #include <stdlib.h>
65
66 #include <vips/vips.h>
67 #include <vips/internal.h>
68 #include "pfreqfilt.h"
69
70 #ifdef HAVE_FFTW
71
72 #include <fftw3.h>
73
74 typedef struct _VipsInvfft {
75 VipsFreqfilt parent_instance;
76
77 gboolean real;
78
79 } VipsInvfft;
80
81 typedef VipsFreqfiltClass VipsInvfftClass;
82
83 G_DEFINE_TYPE( VipsInvfft, vips_invfft, VIPS_TYPE_FREQFILT );
84
85 /* Complex to complex inverse transform.
86 */
87 static int
cinvfft1(VipsObject * object,VipsImage * in,VipsImage ** out)88 cinvfft1( VipsObject *object, VipsImage *in, VipsImage **out )
89 {
90 VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
91 VipsInvfft *invfft = (VipsInvfft *) object;
92 VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( invfft );
93
94 fftw_plan plan;
95 double *planner_scratch;
96
97 if( vips_check_mono( class->nickname, in ) ||
98 vips_check_uncoded( class->nickname, in ) )
99 return( -1 );
100
101 /* Convert input to a complex double membuffer.
102 */
103 *out = vips_image_new_memory();
104 if( vips_cast_dpcomplex( in, &t[0], NULL ) ||
105 vips_image_write( t[0], *out ) )
106 return( -1 );
107
108 /* Make the plan for the transform. Yes, they really do use nx for
109 * height and ny for width.
110 */
111 if( !(planner_scratch = VIPS_ARRAY( invfft,
112 VIPS_IMAGE_N_PELS( in ) * 2, double )) )
113 return( -1 );
114 if( !(plan = fftw_plan_dft_2d( in->Ysize, in->Xsize,
115 (fftw_complex *) planner_scratch,
116 (fftw_complex *) planner_scratch,
117 FFTW_BACKWARD,
118 0 )) ) {
119 vips_error( class->nickname,
120 "%s", _( "unable to create transform plan" ) );
121 return( -1 );
122 }
123
124 fftw_execute_dft( plan,
125 (fftw_complex *) (*out)->data, (fftw_complex *) (*out)->data );
126
127 fftw_destroy_plan( plan );
128
129 (*out)->Type = VIPS_INTERPRETATION_B_W;
130
131 return( 0 );
132 }
133
134 /* Complex to real inverse transform.
135 */
136 static int
rinvfft1(VipsObject * object,VipsImage * in,VipsImage ** out)137 rinvfft1( VipsObject *object, VipsImage *in, VipsImage **out )
138 {
139 VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
140 VipsInvfft *invfft = (VipsInvfft *) object;
141 VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( invfft );
142 const int half_width = in->Xsize / 2 + 1;
143
144 double *half_complex;
145 double *planner_scratch;
146 fftw_plan plan;
147 int x, y;
148 double *q, *p;
149
150 /* Convert input to a complex double membuffer.
151 */
152 t[1] = vips_image_new_memory();
153 if( vips_cast_dpcomplex( in, &t[0], NULL ) ||
154 vips_image_write( t[0], t[1] ) )
155 return( -1 );
156
157 /* Build half-complex image.
158 */
159 if( !(half_complex = VIPS_ARRAY( invfft,
160 t[1]->Ysize * half_width * 2, double )) )
161 return( -1 );
162 q = half_complex;
163 for( y = 0; y < t[1]->Ysize; y++ ) {
164 p = ((double *) t[1]->data) + (guint64) y * t[1]->Xsize * 2;
165
166 for( x = 0; x < half_width; x++ ) {
167 q[0] = p[0];
168 q[1] = p[1];
169 p += 2;
170 q += 2;
171 }
172 }
173
174 /* Make mem buffer real image for output.
175 */
176 *out = vips_image_new_memory();
177 if( vips_image_pipelinev( *out, VIPS_DEMAND_STYLE_ANY, t[1], NULL ) )
178 return( -1 );
179 (*out)->BandFmt = VIPS_FORMAT_DOUBLE;
180 (*out)->Type = VIPS_INTERPRETATION_B_W;
181 if( vips_image_write_prepare( *out ) )
182 return( -1 );
183
184 /* Make the plan for the transform. Yes, they really do use nx for
185 * height and ny for width.
186 */
187 if( !(planner_scratch = VIPS_ARRAY( invfft,
188 t[1]->Ysize * half_width * 2, double )) )
189 return( -1 );
190 if( !(plan = fftw_plan_dft_c2r_2d( t[1]->Ysize, t[1]->Xsize,
191 (fftw_complex *) planner_scratch, (double *) (*out)->data,
192 0 )) ) {
193 vips_error( class->nickname,
194 "%s", _( "unable to create transform plan" ) );
195 return( -1 );
196 }
197
198 fftw_execute_dft_c2r( plan,
199 (fftw_complex *) half_complex, (double *) (*out)->data );
200
201 fftw_destroy_plan( plan );
202
203 return( 0 );
204 }
205
206 static int
vips_invfft_build(VipsObject * object)207 vips_invfft_build( VipsObject *object )
208 {
209 VipsFreqfilt *freqfilt = VIPS_FREQFILT( object );
210 VipsInvfft *invfft = (VipsInvfft *) object;
211 VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
212
213 VipsImage *in;
214
215 if( VIPS_OBJECT_CLASS( vips_invfft_parent_class )->
216 build( object ) )
217 return( -1 );
218
219 in = freqfilt->in;
220
221 if( vips_image_decode( in, &t[0] ) )
222 return( -1 );
223 in = t[0];
224
225 if( invfft->real ) {
226 if( vips__fftproc( VIPS_OBJECT( invfft ),
227 in, &t[1], rinvfft1 ) )
228 return( -1 );
229 }
230 else {
231 if( vips__fftproc( VIPS_OBJECT( invfft ),
232 in, &t[1], cinvfft1 ) )
233 return( -1 );
234 }
235
236 if( vips_image_write( t[1], freqfilt->out ) )
237 return( -1 );
238
239 return( 0 );
240 }
241
242 static void
vips_invfft_class_init(VipsInvfftClass * class)243 vips_invfft_class_init( VipsInvfftClass *class )
244 {
245 GObjectClass *gobject_class = G_OBJECT_CLASS( class );
246 VipsObjectClass *vobject_class = VIPS_OBJECT_CLASS( class );
247
248 gobject_class->set_property = vips_object_set_property;
249 gobject_class->get_property = vips_object_get_property;
250
251 vobject_class->nickname = "invfft";
252 vobject_class->description = _( "inverse FFT" );
253 vobject_class->build = vips_invfft_build;
254
255 VIPS_ARG_BOOL( class, "real", 4,
256 _( "Real" ),
257 _( "Output only the real part of the transform" ),
258 VIPS_ARGUMENT_OPTIONAL_INPUT,
259 G_STRUCT_OFFSET( VipsInvfft, real ),
260 FALSE );
261
262 }
263
264 static void
vips_invfft_init(VipsInvfft * invfft)265 vips_invfft_init( VipsInvfft *invfft )
266 {
267 }
268
269 #endif /*HAVE_FFTW*/
270
271 /**
272 * vips_invfft: (method)
273 * @in: input image
274 * @out: (out): output image
275 * @...: %NULL-terminated list of optional named arguments
276 *
277 * Optional arguments:
278 *
279 * * @real: only output the real part
280 *
281 * Transform an image from Fourier space to real space. The result is complex.
282 * If you are OK with a real result, set @real, it's quicker.
283 *
284 * VIPS uses the fftw Fourier Transform library. If this library was not
285 * available when VIPS was configured, these functions will fail.
286 *
287 * See also: vips_fwfft().
288 *
289 * Returns: 0 on success, -1 on error.
290 */
291 int
vips_invfft(VipsImage * in,VipsImage ** out,...)292 vips_invfft( VipsImage *in, VipsImage **out, ... )
293 {
294 va_list ap;
295 int result;
296
297 va_start( ap, out );
298 result = vips_call_split( "invfft", ap, in, out );
299 va_end( ap );
300
301 return( result );
302 }
303
304