im_fwfft() is a class

ChangeLog

@@ -32,6 +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
 
 20/11/13	started 7.36.5
 - better cache sizing in unbuffered sequential mode

libvips/freq_filt/Makefile.am

@@ -1,6 +1,9 @@
 noinst_LTLIBRARIES = libfreq_filt.la
 
 libfreq_filt_la_SOURCES = \
+	freqfilt.c \
+	pfreqfilt.h \
+	fwfft.c \
 	im_phasecor_fft.c \
 	freq_dispatch.c \
 	im_disp_ps.c \

libvips/freq_filt/freq_dispatch.c

@@ -42,17 +42,6 @@
 
 #include <vips/vips.h>
 
-/**
- * SECTION: freq_filt
- * @short_description: fourier transforms and frequency-domin filters
- * @stability: Stable
- * @see_also: <link linkend="libvips-image">image</link>
- * @include: vips/vips.h
- *
- * To and from Fourier space, filter in Fourier space, convert Fourier-space
- * images to a displayable form.
- */
-
 /* One image in, one out.
  */
 static im_arg_desc one_in_one_out[] = {

libvips/freq_filt/freqfilt.c

@@ -0,0 +1,166 @@
+/* base class for all Fourier stuff 
+ *
+ * properties:
+ * 	- single output image
+ */
+
+/*
+
+    Copyright (C) 1991-2005 The National Gallery
+
+    This library 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.1 of the License, or (at your option) any later version.
+
+    This library 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 library; 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
+
+ */
+
+/*
+#define DEBUG
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif /*HAVE_CONFIG_H*/
+#include <vips/intl.h>
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include <vips/vips.h>
+#include <vips/internal.h>
+
+#include "pfreqfilt.h"
+
+/**
+ * SECTION: freq_filt
+ * @short_description: fourier transforms and frequency-domin filters
+ * @stability: Stable
+ * @see_also: <link linkend="libvips-image">image</link>
+ * @include: vips/vips.h
+ *
+ * To and from Fourier space, filter in Fourier space, convert Fourier-space
+ * images to a displayable form.
+ */
+
+G_DEFINE_ABSTRACT_TYPE( VipsFreqfilt, vips_freqfilt, VIPS_TYPE_OPERATION );
+
+static int
+vips_freqfilt_build( VipsObject *object )
+{
+	VipsFreqfilt *freqfilt = VIPS_FREQFILT( object );
+
+#ifdef DEBUG
+	printf( "vips_freqfilt_build: " );
+	vips_object_print_name( object );
+	printf( "\n" );
+#endif /*DEBUG*/
+
+	g_object_set( freqfilt, "out", vips_image_new(), NULL ); 
+
+	if( VIPS_OBJECT_CLASS( vips_freqfilt_parent_class )->build( object ) )
+		return( -1 );
+
+	return( 0 );
+}
+
+static void
+vips_freqfilt_class_init( VipsFreqfiltClass *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 = "freqfilt";
+	vobject_class->description = _( "frequency-domain filter operations" );
+	vobject_class->build = vips_freqfilt_build;
+
+	VIPS_ARG_IMAGE( class, "out", 1, 
+		_( "Output" ), 
+		_( "Output image" ),
+		VIPS_ARGUMENT_REQUIRED_OUTPUT, 
+		G_STRUCT_OFFSET( VipsFreqfilt, out ) );
+}
+
+static void
+vips_freqfilt_init( VipsFreqfilt *freqfilt )
+{
+}
+
+/* Transform an n-band image with a 1-band processing function.
+ *
+ * Memory strategy: we need memory buffers for the input and the output of
+ * fftw. In some modes fftw generates only half the output and we construct
+ * the rest.
+ *
+ * input pipeline -> 
+ *   bandsplit -> 
+ *     full memory image, freed when im_*fft*() exits -> 
+ *       fftw ->
+ *         half memory image, freed when im_*fft*() exits ->
+ *           full memory image, freed when @out is freed ->
+ *             partial bandjoin ->
+ *               output pipeline
+ *
+ * im__fftproc() needs to just call im__fftproc_fn directly for 1 band images,
+ * so we can't cache the output in this fn.
+ */
+int 
+vips__fftproc( VipsObject *context, 
+	VipsImage *in, VipsImage **out, VipsFftProcessFn fn )
+{
+	VipsImage **bands = (VipsImage **) 
+		vips_object_local_array( context, in->Bands );
+	VipsImage **fft = (VipsImage **) 
+		vips_object_local_array( context, in->Bands );
+
+	int b;
+
+	if( in->Bands == 1 ) 
+		return( fn( context, in, out ) );
+
+	for( b = 0; b < in->Bands; b++ ) 
+		if( vips_extract_band( in, &bands[b], b, NULL ) ||
+			fn( context, bands[b], &fft[b] ) )
+			return( -1 );
+
+	if( vips_bandjoin( fft, out, in->Bands, NULL ) )
+		return( -1 );
+
+	return( 0 );
+}
+
+/* Called from iofuncs to init all operations in this dir. Use a plugin system
+ * instead?
+ */
+void
+vips_freqfilt_operation_init( void )
+{
+#ifdef HAVE_FFTW
+	extern GType vips_fwfft_get_type( void ); 
+#endif /*HAVE_FFTW*/
+
+#ifdef HAVE_FFTW
+	vips_fwfft_get_type(); 
+#endif /*HAVE_FFTW*/
+}
+

libvips/freq_filt/fwfft.c

@@ -0,0 +1,375 @@
+/* forward FFT
+ *
+ * Author: Nicos Dessipris 
+ * Written on: 12/04/1990
+ * Modified on : 09/05/1990	to cope with float input
+ * Modified on : 08/03/1991 	history removed
+ * Modified on : 03/04/1991	to cope with any input
+ *
+ * 28/6/95 JC
+ *	- rewritten to use im_clip2f() rather than own code
+ * 	- memory leaks fixed
+ * 10/9/98 JC
+ *	- frees memory more quickly
+ * 2/4/02 JC
+ *	- fftw code added
+ * 13/7/02 JC
+ *	- output Type set to IM_TYPE_FOURIER to help nip
+ * 27/2/03 JC
+ *	- exploits real_to_complex() path in libfftw for real input (thanks
+ *	  Matt) for a 2x speed-up
+ * 17/11/03 JC
+ *	- fix a segv for wider than high images in the real_to_complex() path
+ *	  (thanks Andrey)
+ *	- fixes to real_to_complex() path to give the correct result for
+ *	  non-square images, including odd widths and heights
+ * 3/11/04
+ *	- added fftw3 support
+ * 7/2/10
+ * 	- cleanups
+ * 	- gtkdoc
+ * 25/3/10
+ * 	- have a "t" image linked to out to keep the image alive for longer
+ * 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 _VipsFwfft {
+	VipsFreqfilt parent_instance;
+
+	VipsImage *in;
+} VipsFwfft;
+
+typedef VipsFreqfiltClass VipsFwfftClass;
+
+G_DEFINE_TYPE( VipsFwfft, vips_fwfft, VIPS_TYPE_FREQFILT );
+
+/* Real to complex forward transform.
+ */
+static int 
+rfwfft1( VipsFwfft *fwfft, VipsImage *in, VipsImage **out )
+{
+	VipsImage **t = (VipsImage **) 
+		vips_object_local_array( VIPS_OBJECT( fwfft ), 4 );
+	VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( fwfft );
+	const guint64 size = VIPS_IMAGE_N_PELS( in );
+	const int half_width = in->Xsize / 2 + 1;
+
+	double *half_complex;
+	double *planner_scratch;
+
+	fftw_plan plan;
+	double *buf, *q, *p;
+	int x, y;
+
+	if( vips_check_mono( class->nickname, in ) ||
+		vips_check_uncoded( class->nickname, in ) )
+                return( -1 );
+
+	/* Convert input to a real double membuffer.
+	 */
+	t[1] = vips_image_new_buffer();
+	if( vips_cast_double( in, &t[0], NULL ) ||
+		vips_image_write( t[0], t[1] ) ); 
+
+	/* Make the plan for the transform. Yes, they really do use nx for
+	 * height and ny for width. Use a separate scratch buffer for the
+	 * planner, we can't overwrite real->data
+	 */
+	if( !(planner_scratch = VIPS_ARRAY( fwfft, 
+		VIPS_IMAGE_N_PELS( in ), double )) )
+		return( -1 );
+	if( !(half_complex = VIPS_ARRAY( fwfft, 
+		in->Ysize * half_width * 2, double )) )
+		return( -1 );
+	if( !(plan = fftw_plan_dft_r2c_2d( in->Ysize, in->Xsize,
+		planner_scratch, (fftw_complex *) half_complex,
+		0 )) ) {
+                vips_error( class->nickname, 
+			"%s", _( "unable to create transform plan" ) );
+		return( -1 );
+	}
+
+	fftw_execute_dft_r2c( plan,
+		(double *) t[1]->data, (fftw_complex *) half_complex );
+
+	fftw_destroy_plan( plan );
+
+	/* Write to out as another memory buffer. 
+	 */
+	*out = vips_image_new_buffer();
+	if( vips_image_pipelinev( *out, VIPS_DEMAND_STYLE_ANY, in, NULL ) )
+                return( -1 );
+	(*out)->BandFmt = VIPS_FORMAT_DPCOMPLEX;
+	(*out)->Type = VIPS_INTERPRETATION_FOURIER;
+	if( !(buf = VIPS_ARRAY( fwfft, VIPS_IMAGE_N_PELS( *out ), double )) )
+		return( -1 );
+
+	/* Copy and normalise. The right half is the up/down and 
+	 * left/right flip of the left, but conjugated. Do the first 
+	 * row separately, then mirror around the centre row.
+	 */
+	p = half_complex;
+	q = buf;
+
+	for( x = 0; x < half_width; x++ ) {
+		q[0] = p[0] / size;
+		q[1] = p[1] / size;
+		p += 2;
+		q += 2;
+	}
+
+	p = half_complex + ((in->Xsize + 1) / 2 - 1) * 2; 
+
+	for( x = half_width; x < (*out)->Xsize; x++ ) {
+		q[0] = p[0] / size;
+		q[1] = -1.0 * p[1] / size;
+		p -= 2;
+		q += 2;
+	}
+
+	if( vips_image_write_line( *out, 0, (VipsPel *) buf ) )
+		return( -1 );
+
+	for( y = 1; y < (*out)->Ysize; y++ ) {
+		p = half_complex + y * half_width * 2; 
+		q = buf;
+
+		for( x = 0; x < half_width; x++ ) {
+			q[0] = p[0] / size;
+			q[1] = p[1] / size;
+			p += 2;
+			q += 2;
+		}
+
+		/* Good grief. 
+		 */
+		p = half_complex + 2 *
+			(((*out)->Ysize - y + 1) * half_width - 2 + 
+				(in->Xsize & 1));
+
+		for( x = half_width; x < (*out)->Xsize; x++ ) {
+			q[0] = p[0] / size;
+			q[1] = -1.0 * p[1] / size;
+			p -= 2;
+			q += 2;
+		}
+
+		if( vips_image_write_line( *out, y, (VipsPel *) buf ) )
+			return( -1 );
+	}
+
+	return( 0 );
+}
+
+/* Complex to complex forward transform.
+ */
+static int 
+cfwfft1( VipsFwfft *fwfft, VipsImage *in, VipsImage **out )
+{
+	VipsImage **t = (VipsImage **) 
+		vips_object_local_array( VIPS_OBJECT( fwfft ), 4 );
+	VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( fwfft );
+
+	fftw_plan plan;
+	double *planner_scratch;
+	double *buf, *q, *p;
+	int x, y;
+
+	if( vips_check_mono( class->nickname, in ) ||
+		vips_check_uncoded( class->nickname, in ) )
+                return( -1 );
+
+	/* 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] ) ); 
+
+	/* We have to have a separate buffer for the planner to work on.
+	 */
+	if( !(planner_scratch = VIPS_ARRAY( fwfft, 
+		VIPS_IMAGE_N_PELS( in ) * 2, double )) )
+		return( -1 );
+
+	/* Make the plan for the transform.
+	 */
+	if( !(plan = fftw_plan_dft_2d( in->Ysize, in->Xsize,
+		(fftw_complex *) planner_scratch, 
+		(fftw_complex *) planner_scratch,
+		FFTW_FORWARD, 
+		0 )) ) {
+                vips_error( class->nickname, 
+			"%s", _( "unable to create transform plan" ) );
+		return( -1 );
+	}
+
+	fftw_execute_dft( plan,
+		(fftw_complex *) t[1]->data, (fftw_complex *) t[1]->data );
+
+	fftw_destroy_plan( plan );
+
+	/* Write to out as another memory buffer. 
+	 */
+	*out = vips_image_new_buffer();
+	if( vips_image_pipelinev( *out, VIPS_DEMAND_STYLE_ANY, in, NULL ) )
+                return( -1 );
+	(*out)->BandFmt = VIPS_FORMAT_DPCOMPLEX;
+	(*out)->Type = VIPS_INTERPRETATION_FOURIER;
+	if( !(buf = VIPS_ARRAY( fwfft, VIPS_IMAGE_N_PELS( *out ), double )) )
+		return( -1 );
+
+	/* Copy to out, normalise.
+	 */
+	p = (double *) t[1]->data;
+	for( y = 0; y < (*out)->Ysize; y++ ) {
+		guint64 size = VIPS_IMAGE_N_PELS( *out );
+
+		q = buf;
+
+		for( x = 0; x < (*out)->Xsize; x++ ) {
+			q[0] = p[0] / size;
+			q[1] = p[1] / size;
+			p += 2;
+			q += 2;
+		}
+
+		if( vips_image_write_line( *out, y, (VipsPel *) buf ) )
+			return( -1 );
+	}
+
+	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 )
+{
+	VipsFreqfilt *freqfilt = VIPS_FREQFILT( object );
+	VipsFwfft *fwfft = (VipsFwfft *) object;
+	VipsImage **t = (VipsImage **) vips_object_local_array( object, 4 );
+
+	if( VIPS_OBJECT_CLASS( vips_fwfft_parent_class )->
+		build( object ) )
+		return( -1 );
+
+	if( vips__fftproc( VIPS_OBJECT( fwfft ), fwfft->in, &t[0], fwfft1 ) ||
+		vips_image_write( t[0], freqfilt->out ) ) 
+		return( -1 );
+
+	return( 0 );
+}
+
+static void
+vips_fwfft_class_init( VipsFwfftClass *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 = "fwfft";
+	vobject_class->description = _( "forward FFT" );
+	vobject_class->build = vips_fwfft_build;
+
+	VIPS_ARG_IMAGE( class, "in", 0, 
+		_( "in" ), 
+		_( "Input image" ),
+		VIPS_ARGUMENT_REQUIRED_INPUT,
+		G_STRUCT_OFFSET( VipsFwfft, in ) );
+
+}
+
+static void
+vips_fwfft_init( VipsFwfft *fwfft )
+{
+}
+
+#endif /*HAVE_FFTW*/
+
+/**
+ * vips_fwfft:
+ * @in: input image 
+ * @out: output image
+ * @...: %NULL-terminated list of optional named arguments
+ *
+ * Transform an image to Fourier space.
+ *
+ * VIPS uses the fftw Fourier Transform library. If this library was not
+ * available when VIPS was configured, these functions will fail.
+ *
+ * See also: vips_invfft().
+ *
+ * Returns: 0 on success, -1 on error.
+ */
+int
+vips_fwfft( VipsImage *in, VipsImage **out, ... )
+{
+	va_list ap;
+	int result;
+
+	va_start( ap, out );
+	result = vips_call_split( "fwfft", ap, in, out );
+	va_end( ap );
+
+	return( result );
+}
+

libvips/freq_filt/im_fwfft.c

@@ -304,6 +304,7 @@ fwfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
 }
 
 #endif 
+typedef int (*im__fftproc_fn)( VipsImage *, VipsImage *, VipsImage * );
 
 /* Transform an n-band image with a 1-band processing function.
  *

libvips/freq_filt/pfreqfilt.h

@@ -0,0 +1,82 @@
+/* base class for all freqfilt operations
+ */
+
+/*
+
+    Copyright (C) 1991-2005 The National Gallery
+
+    This library 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.1 of the License, or (at your option) any later version.
+
+    This library 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 library; 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
+
+ */
+
+#ifndef VIPS_PFREQFILT_H
+#define VIPS_PFREQFILT_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif /*__cplusplus*/
+
+#include <vips/vector.h>
+
+#define VIPS_TYPE_FREQFILT (vips_freqfilt_get_type())
+#define VIPS_FREQFILT( obj ) \
+	(G_TYPE_CHECK_INSTANCE_CAST( (obj), \
+		VIPS_TYPE_FREQFILT, VipsFreqfilt ))
+#define VIPS_FREQFILT_CLASS( klass ) \
+	(G_TYPE_CHECK_CLASS_CAST( (klass), \
+		VIPS_TYPE_FREQFILT, VipsFreqfiltClass))
+#define VIPS_IS_FREQFILT( obj ) \
+	(G_TYPE_CHECK_INSTANCE_TYPE( (obj), VIPS_TYPE_FREQFILT ))
+#define VIPS_IS_FREQFILT_CLASS( klass ) \
+	(G_TYPE_CHECK_CLASS_TYPE( (klass), VIPS_TYPE_FREQFILT ))
+#define VIPS_FREQFILT_GET_CLASS( obj ) \
+	(G_TYPE_INSTANCE_GET_CLASS( (obj), \
+		VIPS_TYPE_FREQFILT, VipsFreqfiltClass ))
+
+typedef struct _VipsFreqfilt {
+	VipsOperation parent_instance;
+
+	/* All have an output image.
+	 */
+	VipsImage *out;
+
+} VipsFreqfilt;
+
+typedef struct _VipsFreqfiltClass {
+	VipsOperationClass parent_class;
+
+} VipsFreqfiltClass;
+
+GType vips_freqfilt_get_type( void );
+
+typedef int (*VipsFftProcessFn)( VipsObject *, VipsImage *, VipsImage ** );
+
+int vips__fftproc( VipsObject *context, 
+	VipsImage *in, VipsImage **out, VipsFftProcessFn fn );
+
+#ifdef __cplusplus
+}
+#endif /*__cplusplus*/
+
+#endif /*VIPS_PFREQFILT_H*/
+
+

libvips/include/vips/internal.h

@@ -105,6 +105,7 @@ void vips__cache_init( void );
 void vips__type_leak( void );
 
 typedef int (*im__fftproc_fn)( VipsImage *, VipsImage *, VipsImage * );
+int im__fftproc( IMAGE *dummy, IMAGE *in, IMAGE *out, im__fftproc_fn fn );
 
 /* iofuncs
  */
@@ -200,9 +201,6 @@ int im__colour_unary( const char *domain,
 VipsImage **im__insert_base( const char *domain, 
 	VipsImage *in1, VipsImage *in2, VipsImage *out );
 
-int im__fftproc( VipsImage *dummy, 
-	VipsImage *in, VipsImage *out, im__fftproc_fn fn );
-
 int im__find_lroverlap( VipsImage *ref_in, VipsImage *sec_in, VipsImage *out,
         int bandno_in,
         int xref, int yref, int xsec, int ysec,
@@ -287,6 +285,7 @@ void vips_resample_operation_init( void );
 void vips_foreign_operation_init( void );
 void vips_colour_operation_init( void );
 void vips_histogram_operation_init( void );
+void vips_freqfilt_operation_init( void );
 
 guint64 vips__parse_size( const char *size_string );
 

libvips/iofuncs/init.c

@@ -270,6 +270,7 @@ vips__init( const char *argv0 )
 	vips_colour_operation_init();
 	vips_histogram_operation_init();
 	vips_convolution_operation_init();
+	vips_freqfilt_operation_init();
 
 	/* Load up any plugins in the vips libdir. We don't error on failure,
 	 * it's too annoying to have VIPS refuse to start because of a broken