DGtal 1.3.0
Loading...
Searching...
No Matches
testRealFFT.cpp
Go to the documentation of this file.
1
33#ifndef WITH_FFTW3
34 #error You need to have activated FFTW3 (WITH_FFTW3) to include this file.
35#endif
36
37#include <cstddef>
38#include <algorithm>
39#include <complex>
40#include <cmath>
41#include <limits>
42#include <string>
43#include <random>
44
45#include <boost/math/constants/constants.hpp>
46
47#include "DGtalCatch.h"
48#include "ConfigTest.h"
49
50#include "DGtal/math/RealFFT.h"
51#include "DGtal/images/ImageContainerBySTLVector.h"
52#include "DGtal/io/readers/PGMReader.h"
53#include "DGtal/io/writers/PGMWriter.h"
54#include "DGtal/base/BasicFunctors.h"
55
56using namespace DGtal;
57
63template <
64 typename TDomain,
65 typename TValue
66>
68{
69 using FFT = RealFFT< TDomain, TValue >;
70
71 INFO( "Initializing RealFFT." );
72 FFT fft( anImage.domain() );
73
74 INFO( "Copying data from the image." );
75 auto spatial_image = fft.getSpatialImage();
76 std::copy( anImage.cbegin(), anImage.cend(), spatial_image.begin() );
77
78 INFO( "Forward transformation." );
79 fft.forwardFFT( FFTW_ESTIMATE );
80
81 INFO( "Checking modification of the input image." );
82 REQUIRE( ! std::equal( anImage.cbegin(), anImage.cend(), spatial_image.cbegin() ) );
83
84 INFO( "Backward transformation." );
85 fft.backwardFFT( FFTW_ESTIMATE );
86
87 INFO( "Comparing result with original image." );
88 const auto eps = 100 * std::numeric_limits<TValue>::epsilon() * std::log( anImage.domain().size() );
89 CAPTURE( eps );
90
91 for ( auto it = spatial_image.cbegin(); it != spatial_image.cend(); ++it )
92 {
93 if ( std::abs( *it - anImage( it.getPoint() ) ) > eps * std::max( *it, TValue(1) ) )
94 FAIL( "Approximation failed: " << *it << " - " << anImage( it.getPoint() ) << " = " << (*it - anImage( it.getPoint() ) ) );
95 }
96}
97
104template <typename TIterator>
105TIterator findMaxNorm( TIterator it, TIterator const & it_end )
106{
107 auto norm_max = std::norm(*it);
108 auto it_max = it;
109
110 for ( ; it != it_end; ++it )
111 if ( std::norm(*it) > norm_max )
112 {
113 norm_max = std::norm(*it);
114 it_max = it;
115 }
116
117 return it_max;
118}
119
125template <
126 typename TDomain,
127 typename TValue
128>
130{
131 typedef RealFFT< TDomain, TValue > FFT; // "typedef" instead of "using" because of g++ 4.7.4 bug.
132 using RealPoint = typename FFT::RealPoint;
133
134 const TValue pi = boost::math::constants::pi<TValue>();
135 const TValue freq = 5;
136 const TValue phase = pi/4;
137
138 INFO( "Checking image size." );
139 REQUIRE( anImage.extent()[ TDomain::dimension-1 ] >= 2*std::abs(freq) );
140
141 INFO( "Initializing RealFFT." );
142 FFT fft( anImage.domain(), RealPoint::zero, RealPoint::diagonal(1) );
143
144 INFO( "Initializing spatial data." );
145 auto spatial_image = fft.getSpatialImage();
146 for ( auto it = spatial_image.begin(); it != spatial_image.end(); ++it )
147 {
148 const auto pt = fft.calcScaledSpatialCoords( it.getPoint() );
149 REQUIRE( fft.calcNativeSpatialCoords( pt ) == it.getPoint() );
150
151 *it = std::cos( 2*pi * freq * pt[ TDomain::dimension - 1 ] + phase );
152 }
153
154 INFO( "Forward transformation." );
155 fft.forwardFFT( FFTW_ESTIMATE );
156
157 INFO( "Finding maximal frequency..." );
158 const auto freq_image = fft.getFreqImage();
159 const auto it_max = findMaxNorm( freq_image.cbegin(), freq_image.cend() );
160 const auto pt_max = it_max.getPoint();
161 INFO( "\tat " << pt_max << " with value " << *it_max );
162
163 INFO( "Checks maximal frequency on unit domain." );
164 {
165 auto freq_pt = fft.calcScaledFreqCoords( it_max.getPoint() );
166 auto freq_val = fft.calcScaledFreqValue( freq_pt, *it_max );
167
168 bool applyConj;
169 REQUIRE( fft.calcNativeFreqCoords( freq_pt, applyConj ) == it_max.getPoint() );
170 REQUIRE( ! applyConj );
171 REQUIRE( std::norm( fft.calcNativeFreqValue( freq_pt, freq_val ) - *it_max ) == Approx(0.).scale(std::norm(*it_max)) );
172 REQUIRE( std::norm( fft.getScaledFreqValue( freq_pt ) - freq_val ) == Approx(0.).scale(std::norm(freq_val)) );
173
174 if ( freq_pt[ TDomain::dimension-1 ] * freq < 0 )
175 {
176 freq_pt = -freq_pt;
177 freq_val = std::conj( freq_val );
178 }
179
180 REQUIRE( ( freq_pt - RealPoint::base( TDomain::dimension-1, freq ) ).norm() == Approx( 0 ).scale(freq_pt.norm()) );
181 REQUIRE( ( std::fmod( std::fmod( std::arg( freq_val ) - phase, 2*pi ) + 3*pi, 2*pi ) - pi ) == Approx( 0 ).scale(pi) );
182
183 }
184
185 INFO( "Checks maximal frequency on translated unit domain." );
186 {
187 fft.setScaledSpatialLowerBound( RealPoint::diagonal( 1. / (4*freq) ) );
188
189 auto freq_pt = fft.calcScaledFreqCoords( it_max.getPoint() );
190 auto freq_val = fft.calcScaledFreqValue( freq_pt, *it_max );
191
192 bool applyConj;
193 REQUIRE( fft.calcNativeFreqCoords( freq_pt, applyConj ) == it_max.getPoint() );
194 REQUIRE( ! applyConj );
195 REQUIRE( std::norm( fft.calcNativeFreqValue( freq_pt, freq_val ) - *it_max ) == Approx(0.).scale(std::norm(*it_max)) );
196 REQUIRE( std::norm( fft.getScaledFreqValue( freq_pt ) - freq_val ) == Approx(0.).scale(std::norm(freq_val)) );
197
198 if ( freq_pt[ TDomain::dimension-1 ] * freq < 0 )
199 {
200 freq_pt = -freq_pt;
201 freq_val = std::conj( freq_val );
202 }
203
204 REQUIRE( ( freq_pt - RealPoint::base( TDomain::dimension-1, freq ) ).norm() == Approx( 0 ).scale(freq_pt.norm()) );
205 REQUIRE( ( std::fmod( std::fmod( std::arg( freq_val ) - phase + pi/2, 2*pi) + 3*pi, 2*pi ) - pi ) == Approx( 0 ).scale(pi) );
206 }
207
208 INFO( "Checks maximal frequency on translated initial domain." );
209 {
210 const RealPoint shift = RealPoint::diagonal( 3. );
211
212 fft.setScaledSpatialExtent( anImage.extent() );
213 fft.setScaledSpatialLowerBound( shift );
214
215 auto freq_pt = fft.calcScaledFreqCoords( it_max.getPoint() );
216 auto freq_val = fft.calcScaledFreqValue( freq_pt, *it_max );
217
218 bool applyConj;
219 REQUIRE( fft.calcNativeFreqCoords( freq_pt, applyConj ) == it_max.getPoint() );
220 REQUIRE( ! applyConj );
221 REQUIRE( std::norm( fft.calcNativeFreqValue( freq_pt, freq_val ) - *it_max ) == Approx(0.).scale(std::norm(*it_max)) );
222 REQUIRE( std::norm( fft.getScaledFreqValue( freq_pt ) - freq_val ) == Approx(0.).scale(std::norm(freq_val)) );
223
224 if ( freq_pt[ TDomain::dimension-1 ] * freq < 0 )
225 {
226 freq_pt = -freq_pt;
227 freq_val = std::conj( freq_val );
228 }
229
230 const auto scaled_factor = freq/anImage.extent()[ TDomain::dimension-1 ];
231 REQUIRE( ( freq_pt - RealPoint::base( TDomain::dimension-1, scaled_factor ) ).norm() == Approx( 0 ).scale(freq_pt.norm()) );
232 REQUIRE( ( std::fmod( std::fmod( std::arg( freq_val ) - phase + 2*pi*scaled_factor*shift[TDomain::dimension-1], 2*pi ) + 3*pi, 2*pi ) - pi ) == Approx( 0 ).scale(pi) );
233 }
234}
235
241template <
242 typename TDomain,
243 typename TValue
244>
246{
247 using FFT = RealFFT< TDomain, TValue >;
248 using Point = typename TDomain::Point;
249
250 const Point shift = anImage.extent() / 3;
251 const auto domain = anImage.domain();
252 const TDomain shifted_domain = TDomain( domain.lowerBound() + shift, domain.upperBound() + shift );
253
254 INFO( "Initializing RealFFT." );
255 FFT fft( domain );
256 //FFT shifted_fft( domain, domain.lowerBound() + shift, anImage.extent() );
257 FFT shifted_fft( shifted_domain );
258
259 INFO( "Pre-creating plan." );
260 fft.createPlan( FFTW_MEASURE, FFTW_FORWARD );
261
262 INFO( "Copying data from the image." );
263 auto spatial_image = fft.getSpatialImage();
264 std::copy( anImage.cbegin(), anImage.cend(), spatial_image.begin() );
265
266 auto shifted_spatial_image = shifted_fft.getSpatialImage();
267 const auto spatial_extent = shifted_fft.getSpatialExtent();
268 for ( auto it = shifted_spatial_image.begin(); it != shifted_spatial_image.end(); ++it )
269 {
270 // Calculating shifted coordinates with periodicity.
271 Point pt = it.getPoint();
272 for ( typename Point::Dimension i = 0; i < Point::dimension; ++i )
273 if ( pt[ i ] > domain.upperBound()[ i ] )
274 pt[ i ] -= anImage.extent()[ i ];
275
276 *it = anImage( pt );
277 }
278
279 INFO( "Forward transformation (forcing plan re-use)." );
280 fft.forwardFFT( FFTW_MEASURE | FFTW_WISDOM_ONLY );
281 shifted_fft.forwardFFT( FFTW_MEASURE | FFTW_WISDOM_ONLY );
282
283 INFO( "Comparing results." );
284 auto freq_image = fft.getFreqImage();
285 auto shifted_freq_image = shifted_fft.getFreqImage();
286 const TValue eps = 100 * std::numeric_limits<TValue>::epsilon();
287
288 for ( auto it = freq_image.cbegin(), shifted_it = shifted_freq_image.cbegin(); it != freq_image.cend(); ++it, ++shifted_it )
289 {
290 if ( std::norm(
291 fft.calcScaledFreqValue( fft.calcScaledFreqCoords( it.getPoint() ), *it )
292 - shifted_fft.calcScaledFreqValue( shifted_fft.calcScaledFreqCoords( shifted_it.getPoint() ), *shifted_it ) )
293 > eps * std::max( std::norm(*it), TValue(1) ) )
294 FAIL( "Approximation failed at point " << it.getPoint()
295 << " between " << *it
296 << " and " << shifted_fft.calcScaledFreqValue( shifted_fft.calcScaledFreqCoords( shifted_it.getPoint() ), *shifted_it )
297 << " (scaled from " << *shifted_it << ")" );
298 }
299
300}
302// Test cases.
303
304#ifdef WITH_FFTW3_FLOAT
305TEST_CASE( "Checking RealFFT on a 2D image in float precision.", "[2D][float]" )
306{
307 constexpr typename DGtal::Dimension N = 2;
308 const std::string file_name = testPath + "/samples/church-small.pgm";
309
310 using real = float;
311 using Space = SpaceND<N>;
314
315 INFO( "Importing image " );
316 const auto image = PGMReader< Image, functors::Cast<real> >::importPGM( file_name );
317
318 INFO( "Testing forward and backward FFT." );
319 testForwardBackwardFFT( image );
320
321 INFO( "Testing spatial domain scaling." );
322 testFFTScaling( image );
323
324 INFO( "Testing FFT on translated image." );
325 cmpTranslatedFFT( image );
326}
327#endif
328
329#ifdef WITH_FFTW3_DOUBLE
330TEST_CASE( "Checking RealFFT on a 2D image in double precision.", "[2D][double]" )
331{
332 constexpr typename DGtal::Dimension N = 2;
333 const std::string file_name = testPath + "/samples/church-small.pgm";
334
335 using real = double;
336 using Space = SpaceND<N>;
339
340 INFO( "Importing image " );
341 const auto image = PGMReader< Image, functors::Cast<real> >::importPGM( file_name );
342
343 INFO( "Testing forward and backward FFT." );
344 testForwardBackwardFFT( image );
345
346 INFO( "Testing spatial domain scaling." );
347 testFFTScaling( image );
348
349 INFO( "Testing FFT on translated image." );
350 cmpTranslatedFFT( image );
351}
352#endif
353
354#ifdef WITH_FFTW3_LONG
355TEST_CASE( "Checking RealFFT on a 2D image in long double precision.", "[2D][long double]" )
356{
357 using namespace DGtal;
358
359 constexpr typename DGtal::Dimension N = 2;
360 const std::string file_name = testPath + "/samples/church-small.pgm";
361
362 using real = long double;
363 using Space = SpaceND<N>;
366
367 INFO( "Importing image " );
368 const auto image = PGMReader< Image, functors::Cast<real> >::importPGM( file_name );
369
370 INFO( "Testing forward and backward FFT." );
371 testForwardBackwardFFT( image );
372
373 INFO( "Testing spatial domain scaling." );
374 testFFTScaling( image );
375
376 INFO( "Testing FFT on translated image." );
377 cmpTranslatedFFT( image );
378}
379#endif
380
381#ifdef WITH_FFTW3_DOUBLE
382TEST_CASE( "Checking RealFFT on a 3D image in double precision.", "[3D][double]" )
383{
384 constexpr typename DGtal::Dimension N = 3;
385
386 using real = double;
387 using Space = SpaceND<N>;
388 using Point = Space::Point;
391
392 const Domain domain( Point{0, 10, 13}, Point{31, 28, 45} );
393 Image image( domain );
394 auto const extent = image.extent();
395
396 INFO( "Filling the image randomly." );
397 const std::size_t CNT = image.size() / 100;
398 std::random_device rd;
399 std::mt19937 gen( rd() );
400 std::uniform_real_distribution<> dis{};
401 Point pt;
402
403 for ( std::size_t i = 0; i < CNT; ++i )
404 {
405 for ( Dimension d = 0; d < N; ++d )
406 pt[ d ] = domain.lowerBound()[ d ] + std::floor( extent[ d ] * dis(gen) );
407 image.setValue( pt, 1. );
408 }
409
410 INFO( "Testing forward and backward FFT." );
411 testForwardBackwardFFT( image );
412
413 INFO( "Testing spatial domain scaling." );
414 testFFTScaling( image );
415
416 INFO( "Testing FFT on translated image." );
417 cmpTranslatedFFT( image );
418}
419#endif
420
421#ifdef WITH_FFTW3_DOUBLE
422TEST_CASE( "Checking RealFFT on a 4D image in double precision.", "[4D][double]" )
423{
424 constexpr typename DGtal::Dimension N = 4;
425
426 using real = double;
427 using Space = SpaceND<N>;
428 using Point = Space::Point;
431
432 const Domain domain( Point{0, 10, 13, 5}, Point{11, 28, 25, 17} );
433 Image image( domain );
434 auto const extent = image.extent();
435
436 INFO( "Filling the image randomly." );
437 const std::size_t CNT = image.size() / 100;
438 std::random_device rd;
439 std::mt19937 gen( rd() );
440 std::uniform_real_distribution<> dis{};
441 Point pt;
442
443 for ( std::size_t i = 0; i < CNT; ++i )
444 {
445 for ( Dimension d = 0; d < N; ++d )
446 pt[ d ] = domain.lowerBound()[ d ] + std::floor( extent[ d ] * dis(gen) );
447 image.setValue( pt, 1. );
448 }
449
450 INFO( "Testing forward and backward FFT." );
451 testForwardBackwardFFT( image );
452
453 INFO( "Testing spatial domain scaling." );
454 testFFTScaling( image );
455
456 INFO( "Testing FFT on translated image." );
457 cmpTranslatedFFT( image );
458}
459#endif
460
461
const Point & lowerBound() const
const Point & upperBound() const
const Vector & extent() const
const Domain & domain() const
Aim: implements association bewteen points lying in a digital domain and values.
Definition: Image.h:70
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
static Self base(Dimension k, Component val=1)
static Self diagonal(Component val=1)
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:1595
PointVector< dim, Integer > Point
Points in DGtal::SpaceND.
Definition: SpaceND.h:110
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:137
Aim: Import a 2D or 3D using the Netpbm formats (ASCII mode).
Definition: PGMReader.h:98
MyPointD Point
Definition: testClone2.cpp:383
CAPTURE(thicknessHV)
Domain domain
void testForwardBackwardFFT(ImageContainerBySTLVector< TDomain, TValue > const &anImage)
Definition: testRealFFT.cpp:67
TIterator findMaxNorm(TIterator it, TIterator const &it_end)
void cmpTranslatedFFT(ImageContainerBySTLVector< TDomain, TValue > const &anImage)
void testFFTScaling(ImageContainerBySTLVector< TDomain, TValue > const &anImage)
TEST_CASE("Checking RealFFT on a 2D image in float precision.", "[2D][float]")
REQUIRE(domain.isInside(aPoint))