remove the fft fallbacks

libvips had a set of fallback fft routines in case fftw was not
detected at configure time ... but they were terrible

remove them! we are fftw-only now
This commit is contained in:
John Cupitt 2012-01-27 14:48:28 +00:00
parent b081f6a2fe
commit 7330c244a4
9 changed files with 44 additions and 509 deletions

View File

@ -69,6 +69,7 @@
- new Python binding based on gobject-introspection
- only spot options at the end of arg strings
- add vips_cache() as a vips8 operator
- remove the old fft fallbacks
14/1/12 started 7.26.8
- interpolate CLI args were broken (thanks speckins)

3
README
View File

@ -100,9 +100,6 @@ fftw3
If libvips finds this library, it uses it for fourier transforms. It
can also use fftw2, but 3 is faster and more accurate.
If the library is not found, libvips falls back to it's own internal
FFT routines which are slower and less accurate.
lcms2
lcms
If present, im_icc_import(), _export() and _transform() are available

33
TODO
View File

@ -1,41 +1,8 @@
- seems to do a lot of looping?
$ vips --vips-progress im_fwfft Gugg_coloured.jpg x.v
vips fwfft1:1: 2 threads, 972 x 16 tiles, groups of 64 scanlines
vips temp-6: 2 threads, 972 x 16 tiles, groups of 64 scanlines
vips temp-6: done in 0s
vips fwfft1:1: done in 0s
vips fft: 2 threads, 128 x 128 tiles, groups of 256 scanlines
vips fft: done in 0s
vips fwfft1:1: 2 threads, 972 x 16 tiles, groups of 64 scanlines
vips fwfft1:1: done in 0s
vips fft: 2 threads, 128 x 128 tiles, groups of 256 scanlines
vips fft: done in 0s
vips fwfft1:1: 2 threads, 972 x 16 tiles, groups of 64 scanlines
vips fwfft1:1: done in 0s
vips fft: 2 threads, 128 x 128 tiles, groups of 256 scanlines
vips fft: done in 0s
vips im__fftproc: 2 threads, 128 x 128 tiles, groups of 256 scanlines
vips im__fftproc: done in 0s
vips x.v: 2 threads, 128 x 128 tiles, groups of 256 scanlines
vips x.v: done in 0s
that's 8 separate loops over the image
- viewing an image tagged as fourier in nip2 is broken, which operation is
failing?
- add a fft test to the suite
remove the old fft fallback

View File

@ -1,7 +1,6 @@
noinst_LTLIBRARIES = libfreq_filt.la
libfreq_filt_la_SOURCES = \
fft_sp.c \
fmask4th.c \
im_phasecor_fft.c \
fmaskcir.c \

View File

@ -1,209 +0,0 @@
/* Copyright (c) 1982 Michael Landy, Yoav Cohen, and George Sperling
Disclaimer: No guarantees of performance accompany this software,
nor is any responsibility assumed on the part of the authors. All the
software has been tested extensively and every effort has been made to
insure its reliability. */
/* 1991 Modified by N. Dessipris to return a valid code on error */
/* fft -- fast fourier transform, adapted from
** Gonzalez & Wintz p.87.
**
** No division by N is performed.
**
** Timing: rough estimates:
** two-dimensional arrays, (including copying columns
** back and forth):
** for arrays up to 16X16: less than 1 sec.
** for 32X32 arrays: about 3.0 sec.
** for 64X64 arrays: about 9.0 sec.
** for 128X128 arrays: about 31.0 sec.
**
** Calling sequence:
**
** float *rvec,*ivec;
** int loglen,skip;
**
** fft_2d(rvec,ivec,loglen)
** performs a 2-dimensional fft where loglen is the log of the length of
** a side of the array
**
** fft_2dgen(rvec,ivec,logrows,logcols)
** performs a 2-dimensional fft where logrows is the log of the number of
** rows, and logcols is the log of the number of columns
**
** fftn(rvec,ivec,loglen,skip)
** performs a 1-dimensional fft on every skip-th entry
*/
/*
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 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>
/* Only compile this if we're missing the fftw library.
*/
#if !HAVE_FFTW && !HAVE_FFTW3
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vips/vips.h>
static float *Const = NULL, *iConst = NULL;
static int storesize = 0, Constsize = 0;
static int
fftn(rvec,ivec,loglen,nskip)
float *rvec,*ivec;
int loglen,nskip;
{
int n,nv2,nm1,i,j,k,l,le,le1,c,nle;
float *rveci , *rvecj , *iveci , *ivecj ;
float t,wr,wi,tr,ti ;
if(loglen==0)
return(-1);
n=1<<loglen ;
nv2=n >> 1 ; nm1=n-1 ; j=0 ;
if (storesize<nv2) {
if ((0==(Const=(float *)calloc((unsigned)nv2,sizeof(float)))) ||
(0==(iConst=(float *)calloc((unsigned)nv2,sizeof(float))))){
im_errormsg( "Not enough core for fftn" );
return(-1);
}
storesize = nv2;
}
if (Constsize!=nv2) {
Constsize = nv2;
wr = cos(2*IM_PI/n);
wi = -sin(2*IM_PI/n);
Const[0] = 1.;
iConst[0] = 0.;
for (i=1;i<nv2;i++) {
Const[i] = wr*Const[i-1] - wi*iConst[i-1];
iConst[i] = wr*iConst[i-1] + wi*Const[i-1];
}
}
for (i=0;i<nm1;i++) {
if(i<j) {
rveci=rvec+i*nskip ; rvecj=rvec+j*nskip ;
t=(*rvecj) ; *rvecj=(*rveci) ; *rveci=t ;
iveci=ivec+i*nskip ; ivecj=ivec+j*nskip ;
t=(*ivecj) ; *ivecj=(*iveci) ; *iveci=t ;
}
k=nv2 ;
while (k<=j) {
j-=k ; k>>=1 ;
}
j+=k ;
}
le=1 ;
for (l=0;l<loglen;l++) {
le1=le ; le+=le ; c = 0; nle = n/le;
for (j=0;j<le1;j++) {
for (i=j;i<n;i+=le) {
if(i+le1>=n) {
im_warning("index=%d\n",i+le1);
return(-1);
}
rveci=rvec+i*nskip ; rvecj=rvec+(i+le1)*nskip;
iveci=ivec+i*nskip ; ivecj=ivec+(i+le1)*nskip;
if (c==0) {
tr = *rvecj;
ti = *ivecj;
}
else {
tr = *rvecj*Const[c] - *ivecj*iConst[c];
ti = *rvecj*iConst[c] + *ivecj*Const[c];
}
*rvecj = *rveci - tr;
*ivecj = *iveci - ti;
*rveci += tr;
*iveci += ti;
}
c += nle;
}
}
/* Division by n
for(i=0;i<n;i++)
{rvec[i*nskip]/=n ; ivec[i*nskip]/=n ;}
*/
return(0);
}
static void
free_store( void )
{
if( Const ) {
free( Const );
Const = NULL;
}
if( iConst ) {
free( iConst );
iConst = NULL;
}
storesize = 0;
Constsize = 0;
}
int
im__fft_sp( float *rvec, float *ivec, int logrows, int logcols )
{
int i,rows,cols,size;
rows = 1<<logrows;
cols = 1<<logcols;
size = rows * cols;
for (i=0;i<size;i+=cols)
if ( fftn(rvec+i,ivec+i,logcols,1) == -1) {
free_store();
return(-1);
}
for (i=0;i<cols;i++)
if ( fftn(rvec+i,ivec+i,logrows,cols) == -1) {
free_store();
return(-1);
}
free_store();
return(0);
}
#endif /*HAVE_FFTW*/

View File

@ -34,6 +34,7 @@
* - 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
*/
/*
@ -82,6 +83,7 @@
#include <vips/internal.h>
#ifdef HAVE_FFTW
/* Call rfftw for a 1 band real image.
*/
static int
@ -269,8 +271,9 @@ fwfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
else
return( rfwfft1( dummy, in, out ) );
}
#else /*!HAVE_FFTW*/
#ifdef HAVE_FFTW3
#elif defined HAVE_FFTW3
/* Real to complex forward transform.
*/
static int
@ -478,96 +481,18 @@ fwfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
else
return( rfwfft1( dummy, in, out ) );
}
#else /*!HAVE_FFTW3*/
/* Transform a 1 band image with vips's built-in fft routine.
*/
#else
static int
fwfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
{
int size = VIPS_IMAGE_N_PELS( in );
int bpx = im_ispoweroftwo( in->Xsize );
int bpy = im_ispoweroftwo( in->Ysize );
float *buf, *q, *p1, *p2;
int x, y;
/* Buffers for real and imaginary parts.
*/
IMAGE *real = im_open_local( dummy, "fwfft1:1", "t" );
IMAGE *imag = im_open_local( dummy, "fwfft1:2", "t" );
/* Temporaries.
*/
IMAGE *t1 = im_open_local( dummy, "fwfft1:3", "p" );
if( !real || !imag || !t1 )
return( -1 );
if( im_pincheck( in ) || im_outcheck( out ) )
return( -1 );
if( in->Coding != IM_CODING_NONE || in->Bands != 1 ||
im_iscomplex( in ) ) {
im_error( "im_fwfft",
"%s", _( "one band non-complex uncoded only" ) );
return( -1 );
}
if( !bpx || !bpy ) {
im_error( "im_fwfft",
"%s", _( "sides must be power of 2" ) );
return( -1 );
}
/* Make sure we have a float input image.
*/
if( im_clip2fmt( in, real, IM_BANDFMT_FLOAT ) )
return( -1 );
/* Make a buffer of 0 floats of the same size for the imaginary part.
*/
if( im_black( t1, in->Xsize, in->Ysize, 1 ) )
return( -1 );
if( im_clip2fmt( t1, imag, IM_BANDFMT_FLOAT ) )
return( -1 );
/* Transform!
*/
if( im__fft_sp( (float *) real->data, (float *) imag->data,
bpx - 1, bpy - 1 ) ) {
im_error( "im_fwfft",
"%s", _( "fft_sp failed" ) );
return( -1 );
}
/* WIO to out.
*/
if( im_cp_desc( out, in ) )
return( -1 );
out->BandFmt = IM_BANDFMT_COMPLEX;
out->Type = IM_TYPE_FOURIER;
if( im_setupout( out ) )
return( -1 );
if( !(buf = (float *) IM_ARRAY( dummy,
IM_IMAGE_SIZEOF_LINE( out ), VipsPel )) )
return( -1 );
/* Gather together real and imag parts. We have to normalise output!
*/
for( p1 = (float *) real->data, p2 = (float *) imag->data,
y = 0; y < out->Ysize; y++ ) {
q = buf;
for( x = 0; x < out->Xsize; x++ ) {
q[0] = *p1++ / size;
q[1] = *p2++ / size;
q += 2;
}
if( im_writeline( y, out, (VipsPel *) buf ) )
return( -1 );
}
return( 0 );
im_error( "im_fwfft",
"%s", _( "vips configured without FFT support" ) );
return( -1 );
}
#endif /*HAVE_FFTW3*/
#endif /*HAVE_FFTW*/
#endif
/* Transform an n-band image with a 1-band processing function.
*/

View File

@ -23,6 +23,7 @@
* - gtkdoc
* 27/1/12
* - better setting of interpretation
* - remove own fft fallback code
*/
/*
@ -71,6 +72,7 @@
#include <vips/internal.h>
#ifdef HAVE_FFTW
/* Call fftw for a 1 band image.
*/
static int
@ -115,8 +117,9 @@ invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
return( 0 );
}
#else /*!HAVE_FFTW*/
#ifdef HAVE_FFTW3
#elif defined HAVE_FFTW3
/* Complex to complex inverse transform.
*/
static int
@ -170,96 +173,18 @@ invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
return( 0 );
}
#else /*!HAVE_FFTW3*/
/* Fall back to VIPS's built-in fft
*/
#else
static int
invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
{
int bpx = im_ispoweroftwo( in->Xsize );
int bpy = im_ispoweroftwo( in->Ysize );
float *buf, *q, *p1, *p2;
int x, y;
/* Buffers for real and imaginary parts.
*/
IMAGE *real = im_open_local( dummy, "invfft1:1", "t" );
IMAGE *imag = im_open_local( dummy, "invfft1:2", "t" );
/* Temps.
*/
IMAGE *t1 = im_open_local( dummy, "invfft1:3", "p" );
IMAGE *t2 = im_open_local( dummy, "invfft1:4", "p" );
if( !real || !imag || !t1 )
return( -1 );
if( im_pincheck( in ) || im_outcheck( out ) )
return( -1 );
if( in->Coding != IM_CODING_NONE ||
in->Bands != 1 || !im_iscomplex( in ) ) {
im_error( "im_invfft",
"%s", _( "one band complex uncoded only" ) );
return( -1 );
}
if( !bpx || !bpy ) {
im_error( "im_invfft",
"%s", _( "sides must be power of 2" ) );
return( -1 );
}
/* Make sure we have a single-precision complex input image.
*/
if( im_clip2fmt( in, t1, IM_BANDFMT_COMPLEX ) )
return( -1 );
/* Extract real and imag parts. We have to complement the imaginary.
*/
if( im_c2real( t1, real ) )
return( -1 );
if( im_c2imag( t1, t2 ) || im_lintra( -1.0, t2, 0.0, imag ) )
return( -1 );
/* Transform!
*/
if( im__fft_sp( (float *) real->data, (float *) imag->data,
bpx - 1, bpy - 1 ) ) {
im_error( "im_invfft",
"%s", _( "fft_sp failed" ) );
return( -1 );
}
/* WIO to out.
*/
if( im_cp_desc( out, in ) )
return( -1 );
out->BandFmt = IM_BANDFMT_COMPLEX;
out->Type = IM_TYPE_B_W;
if( im_setupout( out ) )
return( -1 );
if( !(buf = (float *) IM_ARRAY( dummy,
IM_IMAGE_SIZEOF_LINE( out ), VipsPel )) )
return( -1 );
/* Gather together real and imag parts.
*/
for( p1 = (float *) real->data, p2 = (float *) imag->data,
y = 0; y < out->Ysize; y++ ) {
q = buf;
for( x = 0; x < out->Xsize; x++ ) {
q[0] = *p1++;
q[1] = *p2++;
q += 2;
}
if( im_writeline( y, out, (VipsPel *) buf ) )
return( -1 );
}
return( 0 );
im_error( "im_invfft",
"%s", _( "vips configured without FFT support" ) );
return( -1 );
}
#endif /*HAVE_FFTW3*/
#endif /*HAVE_FFTW*/
#endif
/**
* im_invfft:

View File

@ -9,6 +9,9 @@
* - added fftw3 support
* 7/2/10
* - gtkdoc
* 27/1/12
* - better setting of interpretation
* - remove own fft fallback code
*/
/*
@ -57,6 +60,7 @@
#include <vips/internal.h>
#ifdef HAVE_FFTW
/* Use fftw2.
*/
static int
@ -134,8 +138,9 @@ invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
return( 0 );
}
#else /*!HAVE_FFTW*/
#ifdef HAVE_FFTW3
#elif defined HAVE_FFTW3
/* Complex to real inverse transform.
*/
static int
@ -223,93 +228,18 @@ invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
return( 0 );
}
#else /*!HAVE_FFTW3*/
/* Fall back to vips's built-in fft.
*/
#else
static int
invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
{
int bpx = im_ispoweroftwo( in->Xsize );
int bpy = im_ispoweroftwo( in->Ysize );
float *buf, *q, *p1;
int x, y;
/* Buffers for real and imaginary parts.
*/
IMAGE *real = im_open_local( dummy, "invfft1:1", "t" );
IMAGE *imag = im_open_local( dummy, "invfft1:2", "t" );
/* Temps.
*/
IMAGE *t1 = im_open_local( dummy, "invfft1:3", "p" );
IMAGE *t2 = im_open_local( dummy, "invfft1:4", "p" );
if( !real || !imag || !t1 )
return( -1 );
if( im_pincheck( in ) || im_outcheck( out ) )
return( -1 );
if( in->Coding != IM_CODING_NONE ||
in->Bands != 1 || !im_iscomplex( in ) ) {
im_error( "im_invfft",
"%s", _( "one band complex uncoded only" ) );
return( -1 );
}
if( !bpx || !bpy ) {
im_error( "im_invfft",
"%s", _( "sides must be power of 2" ) );
return( -1 );
}
/* Make sure we have a single-precision complex input image.
*/
if( im_clip2fmt( in, t1, IM_BANDFMT_COMPLEX ) )
return( -1 );
/* Extract real and imag parts. We have to complement the imaginary.
*/
if( im_c2real( t1, real ) )
return( -1 );
if( im_c2imag( t1, t2 ) || im_lintra( -1.0, t2, 0.0, imag ) )
return( -1 );
/* Transform!
*/
if( im__fft_sp( (float *) real->data, (float *) imag->data,
bpx - 1, bpy - 1 ) ) {
im_error( "im_invfft",
"%s", _( "fft_sp failed" ) );
return( -1 );
}
/* WIO to out.
*/
if( im_cp_desc( out, in ) )
return( -1 );
out->BandFmt = IM_BANDFMT_FLOAT;
out->Type = IM_TYPE_B_W;
if( im_setupout( out ) )
return( -1 );
if( !(buf = (float *) IM_ARRAY( dummy,
IM_IMAGE_SIZEOF_LINE( out ), VipsPel )) )
return( -1 );
/* Just write real part.
*/
for( p1 = (float *) real->data, y = 0; y < out->Ysize; y++ ) {
q = buf;
for( x = 0; x < out->Xsize; x++ ) {
q[x] = *p1++;
}
if( im_writeline( y, out, (VipsPel *) buf ) )
return( -1 );
}
return( 0 );
im_error( "im_invfftr",
"%s", _( "vips configured without FFT support" ) );
return( -1 );
}
#endif /*HAVE_FFTW3*/
#endif /*HAVE_FFTW*/
#endif
/**
* im_invfftr:

View File

@ -200,8 +200,8 @@ struct im_col_tab_disp *im_col_make_tables_RGB( VipsImage *im,
struct im_col_display *d );
struct im_col_tab_disp *im_col_display_get_table( struct im_col_display *d );
int im__fft_sp( float *rvec, float *ivec, int logrows, int logcols );
int im__fftproc( VipsImage *dummy, VipsImage *in, VipsImage *out, im__fftproc_fn fn );
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,