clean up fft, reduce memuse

This commit is contained in:
John Cupitt 2012-01-28 11:53:45 +00:00
parent d621cd1f38
commit 4ca2f3ed55
3 changed files with 44 additions and 51 deletions

View File

@ -35,6 +35,8 @@
* 27/1/12
* - better setting of interpretation
* - remove own fft fallback code
* - remove fftw2 path
* - reduce memuse
*/
/*
@ -133,7 +135,8 @@ rfwfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
/* WIO to out.
*/
if( im_cp_desc( out, in ) )
if( im_outcheck( out ) ||
im_cp_desc( out, in ) )
return( -1 );
out->BandFmt = IM_BANDFMT_DPCOMPLEX;
out->Type = IM_TYPE_FOURIER;
@ -215,12 +218,10 @@ cfwfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
vips_check_uncoded( "im_fwfft", in ) )
return( -1 );
/* Double-complex input.
*/
if( !(cmplx = im_open_local( dummy, "fwfft1:1", "t" )) ||
im_pincheck( in ) ||
im_outcheck( out ) )
return( -1 );
if( im_clip2fmt( in, cmplx, IM_BANDFMT_DPCOMPLEX ) )
im_clip2fmt( in, cmplx, IM_BANDFMT_DPCOMPLEX ) )
return( -1 );
/* We have to have a separate buffer for the planner to work on.
@ -241,6 +242,8 @@ cfwfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
return( -1 );
}
if( im_incheck( cmplx ) )
return( -1 );
fftw_execute_dft( plan,
(fftw_complex *) cmplx->data, (fftw_complex *) cmplx->data );
@ -248,7 +251,8 @@ cfwfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
/* WIO to out.
*/
if( im_cp_desc( out, in ) )
if( im_outcheck( out ) ||
im_cp_desc( out, in ) )
return( -1 );
out->BandFmt = IM_BANDFMT_DPCOMPLEX;
out->Type = IM_TYPE_FOURIER;

View File

@ -24,6 +24,8 @@
* 27/1/12
* - better setting of interpretation
* - remove own fft fallback code
* - remove fftw2 path
* - reduce memuse
*/
/*
@ -75,29 +77,25 @@ static int
invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
{
fftw_plan plan;
IMAGE *cmplx;
double *planner_scratch;
IMAGE *cmplx = im_open_local( out, "invfft1:1", "t" );
/* We have to have a separate buffer for the planner to work on.
*/
double *planner_scratch = IM_ARRAY( dummy,
VIPS_IMAGE_N_PELS( in ) * 2, double );
/* Make dp complex image.
*/
if( !cmplx || im_pincheck( in ) || im_poutcheck( out ) )
return( -1 );
if( in->Coding != IM_CODING_NONE || in->Bands != 1 ) {
im_error( "im_invfft",
"%s", _( "one band uncoded only" ) );
if( vips_check_mono( "im_invfft", in ) ||
vips_check_uncoded( "im_invfft", in ) )
return( -1 );
}
if( im_clip2fmt( in, cmplx, IM_BANDFMT_DPCOMPLEX ) )
/* Double-complex input.
*/
if( !(cmplx = im_open_local( dummy, "invfft:1", "t" )) ||
im_clip2fmt( in, cmplx, IM_BANDFMT_DPCOMPLEX ) )
return( -1 );
/* Make the plan for the transform. Yes, they really do use nx for
* height and ny for width.
*/
if( !(planner_scratch = IM_ARRAY( dummy,
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,
@ -108,6 +106,8 @@ invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
return( -1 );
}
if( im_incheck( cmplx ) )
return( -1 );
fftw_execute_dft( plan,
(fftw_complex *) cmplx->data, (fftw_complex *) cmplx->data );
@ -115,8 +115,6 @@ invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
cmplx->Type = IM_TYPE_B_W;
/* Copy to out.
*/
if( im_copy( cmplx, out ) )
return( -1 );

View File

@ -12,6 +12,8 @@
* 27/1/12
* - better setting of interpretation
* - remove own fft fallback code
* - remove fftw2 path
* - reduce memuse
*/
/*
@ -62,41 +64,27 @@
static int
invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
{
IMAGE *cmplx = im_open_local( dummy, "invfft1-1", "t" );
IMAGE *real = im_open_local( out, "invfft1-2", "t" );
const int half_width = in->Xsize / 2 + 1;
/* Transform to halfcomplex here.
*/
double *half_complex = IM_ARRAY( dummy,
in->Ysize * half_width * 2, double );
/* We have to have a separate real buffer for the planner to work on.
*/
double *planner_scratch = IM_ARRAY( dummy,
in->Ysize * half_width * 2, double );
IMAGE *cmplx;
double *half_complex;
IMAGE *real;
double *planner_scratch;
fftw_plan plan;
int x, y;
double *q, *p;
if( !cmplx || !real || !half_complex ||
im_pincheck( in ) ||
im_poutcheck( out ) )
return( -1 );
if( in->Coding != IM_CODING_NONE || in->Bands != 1 ) {
im_error( "im_invfft",
"%s", _( "one band uncoded only" ) );
return( -1 );
}
/* Make dp complex image for input.
/* Double-complex input.
*/
if( im_clip2fmt( in, cmplx, IM_BANDFMT_DPCOMPLEX ) )
if( !(cmplx = im_open_local( dummy, "invfft:1", "t" )) ||
im_clip2fmt( in, cmplx, IM_BANDFMT_DPCOMPLEX ) )
return( -1 );
/* Build half-complex image.
*/
if( !(half_complex = IM_ARRAY( dummy,
in->Ysize * half_width * 2, double )) )
return( -1 );
if( im_incheck( cmplx ) )
return( -1 );
q = half_complex;
@ -113,6 +101,8 @@ invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
/* Make mem buffer real image for output.
*/
if( !(real = im_open_local( out, "invfft1-2", "t" )) )
return( -1 );
if( im_cp_desc( real, in ) )
return( -1 );
real->BandFmt = IM_BANDFMT_DOUBLE;
@ -124,6 +114,9 @@ invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
/* Make the plan for the transform. Yes, they really do use nx for
* height and ny for width.
*/
if( !(planner_scratch = IM_ARRAY( dummy,
in->Ysize * half_width * 2, double )) )
return( -1 );
if( !(plan = fftw_plan_dft_c2r_2d( in->Ysize, in->Xsize,
(fftw_complex *) planner_scratch, (double *) real->data,
0 )) ) {
@ -137,8 +130,6 @@ invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
fftw_destroy_plan( plan );
/* Copy to out.
*/
if( im_copy( real, out ) )
return( -1 );