DGtal  1.5.beta
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 
56 using namespace DGtal;
57 
63 template <
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 
104 template <typename TIterator>
105 TIterator 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 
125 template <
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 
241 template <
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
305 TEST_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
330 TEST_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
355 TEST_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
382 TEST_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
422 TEST_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 Domain & domain() const
const Vector & extent() const
Aim: implements association bewteen points lying in a digital domain and values.
Definition: Image.h:70
float scale
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:136
Aim: Import a 2D or 3D using the Netpbm formats (ASCII mode).
Definition: PGMReader.h:98
int max(int a, int b)
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]")
Image image(domain)
REQUIRE(domain.isInside(aPoint))
PointVector< 3, double > RealPoint