redo im_invfft() as a class

ChangeLog

@@ -32,7 +32,7 @@
 - vipsthumbnail has a --crop option
 - remove video4linux1 code, it was useless on all modern linuxes
 - redone freq filter builders as classes
-- redone im_fwfft() as a class
+- redone im_fwfft(), im_invfft() as classes
 
 20/11/13	started 7.36.5
 - better cache sizing in unbuffered sequential mode

libvips/freq_filt/Makefile.am

@@ -4,6 +4,7 @@ libfreq_filt_la_SOURCES = \
 	freqfilt.c \
 	pfreqfilt.h \
 	fwfft.c \
+	invfft.c \
 	im_phasecor_fft.c \
 	freq_dispatch.c \
 	im_disp_ps.c \

libvips/freq_filt/freqfilt.c

@@ -157,10 +157,12 @@ vips_freqfilt_operation_init( void )
 {
 #ifdef HAVE_FFTW
 	extern GType vips_fwfft_get_type( void ); 
+	extern GType vips_invfft_get_type( void ); 
 #endif /*HAVE_FFTW*/
 
 #ifdef HAVE_FFTW
 	vips_fwfft_get_type(); 
+	vips_invfft_get_type(); 
 #endif /*HAVE_FFTW*/
 }
 

libvips/freq_filt/fwfft.c

@@ -94,10 +94,10 @@ G_DEFINE_TYPE( VipsFwfft, vips_fwfft, VIPS_TYPE_FREQFILT );
 /* Real to complex forward transform.
  */
 static int 
-rfwfft1( VipsFwfft *fwfft, VipsImage *in, VipsImage **out )
+rfwfft1( VipsObject *object, VipsImage *in, VipsImage **out )
 {
-	VipsImage **t = (VipsImage **) 
+	VipsFwfft *fwfft = (VipsFwfft *) object;
-		vips_object_local_array( VIPS_OBJECT( fwfft ), 4 );
+	VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
 	VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( fwfft );
 	const guint64 size = VIPS_IMAGE_N_PELS( in );
 	const int half_width = in->Xsize / 2 + 1;
@@ -212,10 +212,10 @@ rfwfft1( VipsFwfft *fwfft, VipsImage *in, VipsImage **out )
 /* Complex to complex forward transform.
  */
 static int 
-cfwfft1( VipsFwfft *fwfft, VipsImage *in, VipsImage **out )
+cfwfft1( VipsObject *object, VipsImage *in, VipsImage **out )
 {
-	VipsImage **t = (VipsImage **) 
+	VipsFwfft *fwfft = (VipsFwfft *) object;
-		vips_object_local_array( VIPS_OBJECT( fwfft ), 4 );
+	VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
 	VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( fwfft );
 
 	fftw_plan plan;
@@ -288,17 +288,6 @@ cfwfft1( VipsFwfft *fwfft, VipsImage *in, VipsImage **out )
 	return( 0 );
 }
 
-static int 
-fwfft1( VipsObject *object, VipsImage *in, VipsImage **out )
-{
-	VipsFwfft *fwfft = (VipsFwfft *) object;
-
-	if( vips_bandfmt_iscomplex( in->BandFmt ) )
-		return( cfwfft1( fwfft, in, out ) );
-	else
-		return( rfwfft1( fwfft, in, out ) );
-}
-
 static int
 vips_fwfft_build( VipsObject *object )
 {
@@ -310,8 +299,18 @@ vips_fwfft_build( VipsObject *object )
 		build( object ) )
 		return( -1 );
 
-	if( vips__fftproc( VIPS_OBJECT( fwfft ), fwfft->in, &t[0], fwfft1 ) ||
+	if( vips_bandfmt_iscomplex( fwfft->in->BandFmt ) ) {
-		vips_image_write( t[0], freqfilt->out ) ) 
+		if( vips__fftproc( VIPS_OBJECT( fwfft ), fwfft->in, &t[0], 
+			cfwfft1 ) )
+			return( -1 );
+	}
+	else {
+		if( vips__fftproc( VIPS_OBJECT( fwfft ), fwfft->in, &t[0], 
+			rfwfft1 ) )
+			return( -1 );
+	}
+
+	if( vips_image_write( t[0], freqfilt->out ) ) 
 		return( -1 );
 
 	return( 0 );

libvips/freq_filt/invfft.c

@@ -0,0 +1,299 @@
+/* Inverse FFT
+ *
+ * Author: Nicos Dessipris 
+ * Written on: 12/04/1990
+ * Modified on :
+ * 28/6/95 JC
+ *	- rewritten, based on new im_invfft() code
+ * 10/9/98 JC
+ *	- frees memory more quickly
+ * 2/4/02 JC
+ *	- fftw code added
+ * 13/7/02 JC
+ *	- Type reset
+ * 27/2/03 JC
+ *	- tiny speed-up ... save 1 copy on write
+ * 22/1/04 JC
+ *	- oops, fix for segv on wider than high fftw transforms
+ * 3/11/04
+ *	- added fftw3 support
+ * 7/2/10
+ * 	- gtkdoc
+ * 27/1/12
+ * 	- better setting of interpretation
+ * 	- remove own fft fallback code
+ * 	- remove fftw2 path
+ * 	- reduce memuse
+ * 3/1/14
+ * 	- redone as a class
+ */
+
+/*
+
+    This file is part of VIPS.
+    
+    VIPS is free software; you can redistribute it and/or modify
+    it under the terms of the GNU Lesser General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser General Public License for more details.
+
+    You should have received a copy of the GNU Lesser General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+    02110-1301  USA
+
+ */
+
+/*
+
+    These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
+
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif /*HAVE_CONFIG_H*/
+#include <vips/intl.h>
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <vips/vips.h>
+#include "pfreqfilt.h"
+
+#ifdef HAVE_FFTW
+
+#include <fftw3.h>
+
+typedef struct _VipsInvfft {
+	VipsFreqfilt parent_instance;
+
+	VipsImage *in;
+	gboolean real;
+
+} VipsInvfft;
+
+typedef VipsFreqfiltClass VipsInvfftClass;
+
+G_DEFINE_TYPE( VipsInvfft, vips_invfft, VIPS_TYPE_FREQFILT );
+
+/* Complex to complex inverse transform.
+ */
+static int 
+cinvfft1( VipsObject *object, VipsImage *in, VipsImage **out )
+{
+	VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
+	VipsInvfft *invfft = (VipsInvfft *) object;
+	VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( invfft );
+
+	fftw_plan plan;
+	double *planner_scratch;
+
+	if( vips_check_mono( class->nickname, in ) ||
+		vips_check_uncoded( class->nickname, in ) )
+                return( -1 );
+
+	/* Convert input to a complex double membuffer.
+	 */
+	*out = vips_image_new_buffer();
+	if( vips_cast_dpcomplex( in, &t[0], NULL ) ||
+		vips_image_write( t[0], *out ) ); 
+
+	/* Make the plan for the transform. Yes, they really do use nx for
+	 * height and ny for width.
+	 */
+	if( !(planner_scratch = VIPS_ARRAY( invfft, 
+		VIPS_IMAGE_N_PELS( in ) * 2, double )) )
+		return( -1 );
+	if( !(plan = fftw_plan_dft_2d( in->Ysize, in->Xsize,
+		(fftw_complex *) planner_scratch, 
+		(fftw_complex *) planner_scratch,
+		FFTW_BACKWARD, 
+		0 )) ) {
+                vips_error( class->nickname, 
+			"%s", _( "unable to create transform plan" ) );
+		return( -1 );
+	}
+
+	fftw_execute_dft( plan, 
+		(fftw_complex *) (*out)->data, (fftw_complex *) (*out)->data );
+
+	fftw_destroy_plan( plan );
+
+	(*out)->Type = VIPS_INTERPRETATION_B_W;
+
+	return( 0 );
+}
+
+/* Complex to real inverse transform.
+ */
+static int 
+rinvfft1( VipsObject *object, VipsImage *in, VipsImage **out )
+{
+	VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
+	VipsInvfft *invfft = (VipsInvfft *) object;
+	VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( invfft );
+	const int half_width = in->Xsize / 2 + 1;
+
+	double *half_complex;
+	double *planner_scratch;
+	fftw_plan plan;
+	int x, y;
+	double *q, *p;
+
+	/* Convert input to a complex double membuffer.
+	 */
+	t[1] = vips_image_new_buffer();
+	if( vips_cast_dpcomplex( in, &t[0], NULL ) ||
+		vips_image_write( t[0], t[1] ) ); 
+
+	/* Build half-complex image.
+	 */
+	if( !(half_complex = VIPS_ARRAY( invfft, 
+		t[1]->Ysize * half_width * 2, double )) )
+		return( -1 );
+	q = half_complex;
+	for( y = 0; y < t[1]->Ysize; y++ ) {
+		p = ((double *) t[1]->data) + (guint64) y * t[1]->Xsize * 2; 
+
+		for( x = 0; x < half_width; x++ ) {
+			q[0] = p[0];
+			q[1] = p[1];
+			p += 2;
+			q += 2;
+		}
+	}
+
+	/* Make mem buffer real image for output.
+	 */
+	*out = vips_image_new_buffer();
+	if( vips_image_pipelinev( *out, VIPS_DEMAND_STYLE_ANY, t[1], NULL ) )
+                return( -1 );
+	(*out)->BandFmt = VIPS_FORMAT_DOUBLE;
+	(*out)->Type = VIPS_INTERPRETATION_B_W;
+	if( vips_image_write_prepare( *out ) ) 
+		return( -1 ); 
+
+	/* Make the plan for the transform. Yes, they really do use nx for
+	 * height and ny for width.
+	 */
+	if( !(planner_scratch = VIPS_ARRAY( invfft, 
+		t[1]->Ysize * half_width * 2, double )) )
+		return( -1 );
+	if( !(plan = fftw_plan_dft_c2r_2d( t[1]->Ysize, t[1]->Xsize,
+		(fftw_complex *) planner_scratch, (double *) (*out)->data,
+		0 )) ) {
+                vips_error( class->nickname,
+			"%s", _( "unable to create transform plan" ) );
+		return( -1 );
+	}
+
+	fftw_execute_dft_c2r( plan,
+		(fftw_complex *) half_complex, (double *) (*out)->data );
+
+	fftw_destroy_plan( plan );
+
+	return( 0 );
+}
+
+static int
+vips_invfft_build( VipsObject *object )
+{
+	VipsFreqfilt *freqfilt = VIPS_FREQFILT( object );
+	VipsInvfft *invfft = (VipsInvfft *) object;
+	VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
+
+	if( VIPS_OBJECT_CLASS( vips_invfft_parent_class )->
+		build( object ) )
+		return( -1 );
+
+	if( invfft->real ) {
+		if( vips__fftproc( VIPS_OBJECT( invfft ), 
+			invfft->in, &t[0], rinvfft1 ) )
+			return( -1 );
+	}
+	else {
+		if( vips__fftproc( VIPS_OBJECT( invfft ), 
+			invfft->in, &t[0], cinvfft1 ) )
+			return( -1 );
+	}
+	
+	if( vips_image_write( t[0], freqfilt->out ) ) 
+		return( -1 );
+
+	return( 0 );
+}
+
+static void
+vips_invfft_class_init( VipsInvfftClass *class )
+{
+	GObjectClass *gobject_class = G_OBJECT_CLASS( class );
+	VipsObjectClass *vobject_class = VIPS_OBJECT_CLASS( class );
+
+	gobject_class->set_property = vips_object_set_property;
+	gobject_class->get_property = vips_object_get_property;
+
+	vobject_class->nickname = "invfft";
+	vobject_class->description = _( "inverse FFT" );
+	vobject_class->build = vips_invfft_build;
+
+	VIPS_ARG_IMAGE( class, "in", 0, 
+		_( "in" ), 
+		_( "Input image" ),
+		VIPS_ARGUMENT_REQUIRED_INPUT,
+		G_STRUCT_OFFSET( VipsInvfft, in ) );
+
+	VIPS_ARG_BOOL( class, "real", 4, 
+		_( "Real" ), 
+		_( "Ouput only the real part of the transform" ),
+		VIPS_ARGUMENT_OPTIONAL_INPUT,
+		G_STRUCT_OFFSET( VipsInvfft, real ),
+		FALSE );
+
+}
+
+static void
+vips_invfft_init( VipsInvfft *invfft )
+{
+}
+
+#endif /*HAVE_FFTW*/
+
+/**
+ * vips_invfft:
+ * @in: input image 
+ * @out: output image
+ * @...: %NULL-terminated list of optional named arguments
+ *
+ * Optional arguments:
+ *
+ * @real: only output the real part
+ *
+ * Transform an image from Fourier space to real space. The result is complex.
+ * If you are OK with a real result, set @real, it's quicker.
+ *
+ * VIPS uses the fftw Fourier Transform library. If this library was not
+ * available when VIPS was configured, these functions will fail.
+ *
+ * See also: vips_fwfft().
+ *
+ * Returns: 0 on success, -1 on error.
+ */
+int
+vips_invfft( VipsImage *in, VipsImage **out, ... )
+{
+	va_list ap;
+	int result;
+
+	va_start( ap, out );
+	result = vips_call_split( "invfft", ap, in, out );
+	va_end( ap );
+
+	return( result );
+}
+

libvips/include/vips/create.h

@@ -38,8 +38,6 @@
 extern "C" {
 #endif /*__cplusplus*/
 
-void vips_create_operation_init( void );
-
 int vips_black( VipsImage **out, int width, int height, ... )
 	__attribute__((sentinel));
 

libvips/include/vips/freq_filt.h

@@ -38,7 +38,10 @@
 extern "C" {
 #endif /*__cplusplus*/
 
-#include <vips/vips.h>
+int vips_fwfft( VipsImage *in, VipsImage **out, ... )
+	__attribute__((sentinel));
+int vips_invfft( VipsImage *in, VipsImage **out, ... )
+	__attribute__((sentinel));
 
 int im_fwfft( VipsImage *in, VipsImage *out );
 int im_invfft( VipsImage *in, VipsImage *out );

libvips/include/vips/internal.h

@@ -286,6 +286,7 @@ void vips_foreign_operation_init( void );
 void vips_colour_operation_init( void );
 void vips_histogram_operation_init( void );
 void vips_freqfilt_operation_init( void );
+void vips_create_operation_init( void );
 
 guint64 vips__parse_size( const char *size_string );