DGtal  1.5.beta
testRealFFT.cpp File Reference
#include <cstddef>
#include <algorithm>
#include <complex>
#include <cmath>
#include <limits>
#include <string>
#include <random>
#include <boost/math/constants/constants.hpp>
#include "DGtalCatch.h"
#include "ConfigTest.h"
#include "DGtal/math/RealFFT.h"
#include "DGtal/images/ImageContainerBySTLVector.h"
#include "DGtal/io/readers/PGMReader.h"
#include "DGtal/io/writers/PGMWriter.h"
#include "DGtal/base/BasicFunctors.h"
Include dependency graph for testRealFFT.cpp:

Go to the source code of this file.

Functions

template<typename TDomain , typename TValue >
void testForwardBackwardFFT (ImageContainerBySTLVector< TDomain, TValue > const &anImage)
 
template<typename TIterator >
TIterator findMaxNorm (TIterator it, TIterator const &it_end)
 
template<typename TDomain , typename TValue >
void testFFTScaling (ImageContainerBySTLVector< TDomain, TValue > const &anImage)
 
template<typename TDomain , typename TValue >
void cmpTranslatedFFT (ImageContainerBySTLVector< TDomain, TValue > const &anImage)
 
 TEST_CASE ("Checking RealFFT on a 2D image in float precision.", "[2D][float]")
 
 TEST_CASE ("Checking RealFFT on a 2D image in double precision.", "[2D][double]")
 
 TEST_CASE ("Checking RealFFT on a 2D image in long double precision.", "[2D][long double]")
 
 TEST_CASE ("Checking RealFFT on a 3D image in double precision.", "[3D][double]")
 
 TEST_CASE ("Checking RealFFT on a 4D image in double precision.", "[4D][double]")
 

Detailed Description

This program 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 3 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 General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

Author
Roland Denis (rolan.nosp@m.d.de.nosp@m.nis@u.nosp@m.niv-.nosp@m.smb.f.nosp@m.r )
Date
2016/06/02

This file is part of the DGtal library

Definition in file testRealFFT.cpp.

Function Documentation

◆ cmpTranslatedFFT()

template<typename TDomain , typename TValue >
void cmpTranslatedFFT ( ImageContainerBySTLVector< TDomain, TValue > const &  anImage)

Compares the FFT of an image with the FFT of its translation.

Template Parameters
TDomainDomain type.
TValueValue type.
Parameters
anImageThe image from which to take the domain.

Definition at line 245 of file testRealFFT.cpp.

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 }
const Point & lowerBound() const
const Point & upperBound() const
const Domain & domain() const
const Vector & extent() const
int max(int a, int b)
MyPointD Point
Definition: testClone2.cpp:383
Domain domain

References DGtal::ImageContainerBySTLVector< TDomain, TValue >::domain(), domain, DGtal::ImageContainerBySTLVector< TDomain, TValue >::extent(), DGtal::HyperRectDomain< TSpace >::lowerBound(), max(), and DGtal::HyperRectDomain< TSpace >::upperBound().

Referenced by TEST_CASE().

◆ findMaxNorm()

template<typename TIterator >
TIterator findMaxNorm ( TIterator  it,
TIterator const &  it_end 
)

Find element with maximal norm in a range.

Template Parameters
TIteratorIterator type.
Parameters
itIterator to the range begin.
it_endIterator the the range end.
Returns
the iterator to the maximal element.

Definition at line 105 of file testRealFFT.cpp.

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 }

Referenced by testFFTScaling().

◆ TEST_CASE() [1/5]

TEST_CASE ( "Checking RealFFT on a 2D image in double precision."  ,
""  [2D][double] 
)

Definition at line 330 of file testRealFFT.cpp.

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 }
Aim: implements association bewteen points lying in a digital domain and values.
Definition: Image.h:70
DGtal::uint32_t Dimension
Definition: Common.h:136
Aim: Import a 2D or 3D using the Netpbm formats (ASCII mode).
Definition: PGMReader.h:98
void testForwardBackwardFFT(ImageContainerBySTLVector< TDomain, TValue > const &anImage)
Definition: testRealFFT.cpp:67
void cmpTranslatedFFT(ImageContainerBySTLVector< TDomain, TValue > const &anImage)
void testFFTScaling(ImageContainerBySTLVector< TDomain, TValue > const &anImage)

References cmpTranslatedFFT(), testFFTScaling(), and testForwardBackwardFFT().

◆ TEST_CASE() [2/5]

TEST_CASE ( "Checking RealFFT on a 2D image in float precision."  ,
""  [2D][float] 
)

Definition at line 305 of file testRealFFT.cpp.

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 }

References cmpTranslatedFFT(), testFFTScaling(), and testForwardBackwardFFT().

◆ TEST_CASE() [3/5]

TEST_CASE ( "Checking RealFFT on a 2D image in long double precision."  ,
""  [2D][long double] 
)

Definition at line 355 of file testRealFFT.cpp.

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 }
DGtal is the top-level namespace which contains all DGtal functions and types.

References cmpTranslatedFFT(), testFFTScaling(), and testForwardBackwardFFT().

◆ TEST_CASE() [4/5]

TEST_CASE ( "Checking RealFFT on a 3D image in double precision."  ,
""  [3D][double] 
)

Definition at line 382 of file testRealFFT.cpp.

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 }
Image image(domain)

References cmpTranslatedFFT(), domain, image(), DGtal::HyperRectDomain< TSpace >::lowerBound(), testFFTScaling(), and testForwardBackwardFFT().

◆ TEST_CASE() [5/5]

TEST_CASE ( "Checking RealFFT on a 4D image in double precision."  ,
""  [4D][double] 
)

Definition at line 422 of file testRealFFT.cpp.

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 }

References cmpTranslatedFFT(), domain, image(), DGtal::HyperRectDomain< TSpace >::lowerBound(), testFFTScaling(), and testForwardBackwardFFT().

◆ testFFTScaling()

template<typename TDomain , typename TValue >
void testFFTScaling ( ImageContainerBySTLVector< TDomain, TValue > const &  anImage)

Tests spatial domain scaling.

Template Parameters
TDomainDomain type.
TValueValue type.
Parameters
anImageThe image from which to take the domain.

Definition at line 129 of file testRealFFT.cpp.

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 }
float scale
TIterator findMaxNorm(TIterator it, TIterator const &it_end)
REQUIRE(domain.isInside(aPoint))
PointVector< 3, double > RealPoint

References DGtal::ImageContainerBySTLVector< TDomain, TValue >::domain(), DGtal::ImageContainerBySTLVector< TDomain, TValue >::extent(), findMaxNorm(), REQUIRE(), and scale.

Referenced by TEST_CASE().

◆ testForwardBackwardFFT()

template<typename TDomain , typename TValue >
void testForwardBackwardFFT ( ImageContainerBySTLVector< TDomain, TValue > const &  anImage)

Tests of forward FFT followed by a backward FFT.

Template Parameters
TDomainDomain type.
TValueValue type.
Parameters
anImageThe image to transform.

Definition at line 67 of file testRealFFT.cpp.

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 }
CAPTURE(thicknessHV)

References CAPTURE(), DGtal::ImageContainerBySTLVector< TDomain, TValue >::domain(), max(), and REQUIRE().

Referenced by TEST_CASE().