DGtal  1.5.beta
RealFFT.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file
19  * @author Roland Denis (\c roland.denis@univ-smb.fr )
20  * LAboratory of MAthematics - LAMA (CNRS, UMR 5127), University of Savoie, France
21  *
22  * @date 2016/06/01
23  *
24  * Implementation of inline methods defined in RealFFT.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <ostream>
32 #include <stdexcept> // Exceptions
33 #include <new> // std::bad_alloc exception
34 #include <algorithm>
35 #include <cmath>
36 //////////////////////////////////////////////////////////////////////////////
37 
38 ///////////////////////////////////////////////////////////////////////////////
39 // IMPLEMENTATION of inline methods.
40 ///////////////////////////////////////////////////////////////////////////////
41 
42 ///////////////////////////////////////////////////////////////////////////////
43 // ----------------------- Standard services ------------------------------
44 
45 
46 ///////////////////////////////////////////////////////////////////////////////
47 // Interface - public :
48 
49 // Constructor from a domain and scaled lower bound and extent.
50 template <typename TSpace, typename T>
51 inline
52 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
53  RealFFT( Domain const& aDomain,
54  RealPoint const& aLowerBound,
55  RealPoint const& anExtent
56  )
57  : mySpatialDomain{ aDomain }
58  , mySpatialExtent{ mySpatialDomain.upperBound() - mySpatialDomain.lowerBound() + Point::diagonal(1) }
59  , myFreqExtent{ mySpatialExtent / (Point::diagonal(1) + Point::base(0)) + Point::base(0) }
60  , myFreqDomain{ Point::diagonal(0), myFreqExtent - Point::diagonal(1) }
61  , myFullSpatialDomain{ mySpatialDomain.lowerBound(), mySpatialDomain.lowerBound() + myFreqExtent + Point::base(0, myFreqExtent[0]) - Point::diagonal(1) }
62  , myStorage( FFTW::malloc( sizeof(Complex) * myFreqDomain.size() ) )
63 {
64  if ( myStorage == nullptr ) throw std::bad_alloc{};
65 
66  setScaledSpatialLowerBound( aLowerBound );
67  setScaledSpatialExtent( anExtent );
68 }
69 
70 // Constructor from a domain.
71 template <typename TSpace, typename T>
72 inline
73 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
74  RealFFT( Domain const& aDomain )
75  : RealFFT( aDomain
76  , aDomain.lowerBound()
77  , aDomain.upperBound() - aDomain.lowerBound() + Point::diagonal(1)
78  )
79 {
80 }
81 
82 // Destructor
83 template <typename TSpace, typename T>
84 inline
85 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
86  ~RealFFT()
87 {
88  FFTW::free( myStorage );
89 }
90 
91 // Padding when using real datas.
92 template <typename TSpace, typename T>
93 inline
94 std::size_t
95 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
96  getPadding() const noexcept
97 {
98  return 2*myFreqExtent[0] - mySpatialExtent[0];
99 }
100 
101 // Gets mutable spatial storage.
102 template <typename TSpace, typename T>
103 inline
104 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Real*
105 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
106  getSpatialStorage() noexcept
107 {
108  return reinterpret_cast<Real*>(myStorage);
109 }
110 
111 // Gets non-mutable spatial storage.
112 template <typename TSpace, typename T>
113 inline
114 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Real const*
115 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
116  getSpatialStorage() const noexcept
117 {
118 }
119 
120 // Gets mutable spatial image.
121 template <typename TSpace, typename T>
122 inline
123 DGtal::ArrayImageAdapter<
124  typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Real*,
125  typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Domain
126 >
127 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
128  getSpatialImage() noexcept
129 {
130  return { getSpatialStorage(), myFullSpatialDomain, mySpatialDomain };
131 }
132 
133 // Gets non-mutable spatial image.
134 template <typename TSpace, typename T>
135 inline
136 DGtal::ArrayImageAdapter<
137  typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Real const*,
138  typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Domain
139 >
140 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
141  getSpatialImage() const noexcept
142 {
143  return { getSpatialStorage(), myFullSpatialDomain, mySpatialDomain };
144 }
145 
146 // Gets mutable frequential raw storage.
147 template <typename TSpace, typename T>
148 inline
149 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Complex*
150 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
151  getFreqStorage() noexcept
152 {
153  return reinterpret_cast<Complex*>(myStorage);
154 }
155 
156 // Gets non-mutable frequential storage.
157 template <typename TSpace, typename T>
158 inline
159 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Complex const*
160 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
161  getFreqStorage() const noexcept
162 {
163  return reinterpret_cast<Complex const*>(myStorage);
164 }
165 
166 // Gets mutable frequential image.
167 template <typename TSpace, typename T>
168 inline
169 DGtal::ArrayImageAdapter<
170  typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Complex*,
171  typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Domain
172 >
173 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
174  getFreqImage() noexcept
175 {
176  return { getFreqStorage(), getFreqDomain() };
177 }
178 
179 // Get non-mutable frequential image.
180 template <typename TSpace, typename T>
181 inline
182 DGtal::ArrayImageAdapter<
183  typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Complex const*,
184  typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Domain
185 >
186 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
187  getFreqImage() const noexcept
188 {
189  return { getFreqStorage(), getFreqDomain() };
190 }
191 
192 // Get spatial domain.
193 template <typename TSpace, typename T>
194 inline
195 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Domain const&
196 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
197  getSpatialDomain() const noexcept
198 {
199  return mySpatialDomain;
200 }
201 
202 // Get frequential domain.
203 template <typename TSpace, typename T>
204 inline
205 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Domain const&
206 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
207  getFreqDomain() const noexcept
208 {
209  return myFreqDomain;
210 }
211 
212 // Get spatial domain extent.
213 template <typename TSpace, typename T>
214 inline
215 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Point const&
216 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
217  getSpatialExtent() const noexcept
218 {
219  return mySpatialExtent;
220 }
221 
222 // Get frequential domain extent.
223 template <typename TSpace, typename T>
224 inline
225 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Point const&
226 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
227  getFreqExtent() const noexcept
228 {
229  return myFreqExtent;
230 }
231 
232 // Gets the extent of the scaled spatial domain.
233 template <typename TSpace, typename T>
234 inline
235 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::RealPoint
236 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
237  getScaledSpatialExtent() const noexcept
238 {
239  return myScaledSpatialExtent;
240 }
241 
242 // Sets the extent of the scaled spatial domain.
243 template <typename TSpace, typename T>
244 inline
245 void
246 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
247  setScaledSpatialExtent( RealPoint const& anExtent ) noexcept
248 {
249  myScaledSpatialExtent = anExtent;
250 
251  myScaledFreqMag = 1.;
252  for ( Dimension i = 0; i < dimension; ++i )
253  myScaledFreqMag *= mySpatialExtent[ i ] / myScaledSpatialExtent[ i ];
254 }
255 
256 // Gets the lower bound of the scaled spatial domain.
257 template <typename TSpace, typename T>
258 inline
259 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::RealPoint
260 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
261  getScaledSpatialLowerBound() const noexcept
262 {
263  return myScaledSpatialLowerBound;
264 }
265 
266 // Sets the lower bound of the scaled spatial domain.
267 template <typename TSpace, typename T>
268 inline
269 void
270 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
271  setScaledSpatialLowerBound( RealPoint const& aPoint ) noexcept
272 {
273  myScaledSpatialLowerBound = aPoint;
274 }
275 
276 // Converts coordinates from the spatial image into scaled coordinates.
277 template <typename TSpace, typename T>
278 inline
279 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::RealPoint
280 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
281  calcScaledSpatialCoords( Point const& aPoint ) const noexcept
282 {
283  RealPoint scaledCoords;
284  for ( Dimension i = 0; i < dimension; ++i )
285  scaledCoords[ i ] = myScaledSpatialLowerBound[ i ]
286  + ( myScaledSpatialExtent[ i ] / mySpatialExtent[ i ] )
287  * ( aPoint[ i ] - mySpatialDomain.lowerBound()[ i ] );
288 
289  return scaledCoords;
290 }
291 
292 // Converts scaled coordinates from the spatial image into native spatial coords.
293 template <typename TSpace, typename T>
294 inline
295 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Point
296 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
297  calcNativeSpatialCoords( RealPoint const& aScaledPoint ) const noexcept
298 {
299  Point nativeCoords;
300  for ( Dimension i = 0; i < dimension; ++i )
301  nativeCoords[ i ] = std::round(
302  mySpatialDomain.lowerBound()[ i ]
303  + ( mySpatialExtent[ i ] / myScaledSpatialExtent[ i ] )
304  * ( aScaledPoint[ i ] - myScaledSpatialLowerBound[ i ] )
305  );
306 
307  return nativeCoords;
308 }
309 
310 // Converts coordinates from the frequency image into scaled frequencies.
311 template <typename TSpace, typename T>
312 inline
313 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::RealPoint
314 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
315  calcScaledFreqCoords( Point const& aPoint ) const noexcept
316 {
317  RealPoint frequencies;
318  for ( Dimension i = 0; i < dimension; ++i )
319  {
320  if ( aPoint[ i ] <= mySpatialExtent[ i ] / 2 )
321  frequencies[ i ] = aPoint[ i ] / myScaledSpatialExtent[ i ];
322  else
323  frequencies[ i ] = ( aPoint[ i ] - mySpatialExtent[ i ] ) / myScaledSpatialExtent[ i ];
324  }
325 
326  return frequencies;
327 }
328 
329 // From the scaled frequency coordinates, calculates the nearest corresponding
330 // coordinates in the native frequency image.
331 template <typename TSpace, typename T>
332 inline
333 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Point
334 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
335  calcNativeFreqCoords( RealPoint const& aScaledPoint, bool & applyConj ) const noexcept
336 {
337  Point nativeCoords;
338  for ( Dimension i = 0; i < dimension; ++i )
339  {
340  nativeCoords[ i ] = std::round( aScaledPoint[ i ] * myScaledSpatialExtent[ i ] );
341 
342  if ( aScaledPoint[ i ] < 0 )
343  nativeCoords[ i ] += mySpatialExtent[ i ];
344  }
345 
346  // Checking if conjugate is needed.
347  if ( nativeCoords[ 0 ] > mySpatialExtent[ 0 ] / 2 )
348  {
349  for ( Dimension i = 0; i < dimension; ++i )
350  if ( nativeCoords[ i ] > 0 )
351  nativeCoords[ i ] = mySpatialExtent[ i ] - nativeCoords[ i ];
352 
353  applyConj = true;
354  }
355  else
356  applyConj = false;
357 
358  return nativeCoords;
359 }
360 
361 
362 // Converts a complex value from the frequency image to the scaled frequency image.
363 template <typename TSpace, typename T>
364 inline
365 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Complex
366 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
367  calcScaledFreqValue( RealPoint const& aScaledPoint, Complex const& aValue ) const noexcept
368 {
369  Real phase = 0;
370  for ( Dimension i = 0; i < dimension; ++i )
371  phase -= Real( myScaledSpatialLowerBound[ i ] ) * aScaledPoint[ i ];
372 
373  phase *= 2*pi;
374 
375  return aValue * std::polar( myScaledFreqMag, phase );
376 }
377 
378 // From a complex value of the scaled frequency image, calculates the
379 // corresponding value (undoing scaling and rotation) in the native
380 // frequency image.
381 template <typename TSpace, typename T>
382 inline
383 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Complex
384 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
385  calcNativeFreqValue( RealPoint const& aScaledPoint, Complex const& aScaledValue ) const noexcept
386 {
387  Real phase = 0;
388  for ( Dimension i = 0; i < dimension; ++i )
389  phase += Real( myScaledSpatialLowerBound[ i ] ) * aScaledPoint[ i ];
390 
391  phase *= 2*pi;
392 
393  return aScaledValue * std::polar( Real(1)/myScaledFreqMag, phase );
394 }
395 
396 // Creates a transformation plan for the specified transformation direction.
397 template <typename TSpace, typename T>
398 inline
399 void
400 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
401  createPlan( unsigned flags, int way )
402 {
403  // Transform dimensions
404  int n[dimension];
405  for (size_t i = 0; i < dimension; ++i)
406  n[dimension-i-1] = mySpatialExtent[i];
407 
408  // Creates the plan for this transformation
409  typename FFTW::plan p = FFTW::plan_dft( dimension, n, getSpatialStorage(), getFreqStorage(), way, flags );
410 
411  // Destroying plan
412  FFTW::destroy_plan( p );
413 }
414 
415 // Fast Fourier Transformation.
416 template <typename TSpace, typename T>
417 inline
418 void
419 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
420  doFFT( unsigned flags, int way, bool normalized )
421 {
422  typename FFTW::plan p;
423 
424  // Transform dimensions
425  int n[dimension];
426  for (size_t i = 0; i < dimension; ++i)
427  n[dimension-i-1] = mySpatialExtent[i];
428 
429  // Creates the plan for this transformation
430  // Only FFTW_ESTIMATE flag preserves input.
431  if ( flags & FFTW_ESTIMATE )
432  {
433  p = FFTW::plan_dft( dimension, n, getSpatialStorage(), getFreqStorage(), way, FFTW_ESTIMATE );
434  }
435  else
436  {
437  // Strategy to preserve input datas while creating DFT plan:
438  // - Firstly, checks if a plan already exists for this dimensions.
439  p = FFTW::plan_dft( dimension, n, getSpatialStorage(), getFreqStorage(), way, FFTW_WISDOM_ONLY | flags );
440 
441  // - Otherwise, create fake input to create the plan.
442  if ( p == NULL )
443  {
444  void* tmp = FFTW::malloc( sizeof(Complex) * myFreqDomain.size() );
445  if ( tmp == nullptr ) throw std::bad_alloc{};
446  p = FFTW::plan_dft( dimension, n, reinterpret_cast<Real*>(tmp), reinterpret_cast<Complex*>(tmp), way, flags );
447  FFTW::free(tmp);
448  }
449  }
450 
451  // We must have a valid plan now ...
452  if ( p == NULL ) throw std::runtime_error("No valid DFT plan founded.");
453 
454  // Gogogo !
455  FFTW::execute_dft( p, getSpatialStorage(), getFreqStorage(), way );
456 
457  // Destroying plan
458  FFTW::destroy_plan( p );
459 
460  // Normalization
461  if ( way == FFTW_BACKWARD && normalized )
462  {
463  const std::size_t N = getSpatialDomain().size();
464 
465  for ( auto & v : getSpatialImage() )
466  v /= N;
467  }
468 }
469 
470 // In-place forward FFT transformation (spatial -> frequential)
471 template <typename TSpace, typename T>
472 inline
473 void
474 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
475  forwardFFT( unsigned flags )
476 {
477  doFFT( flags, FFTW_FORWARD, false );
478 }
479 
480 // In-place backward FFT transformation (frequential -> spatial)
481 template <typename TSpace, typename T>
482 inline
483 void
484 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
485  backwardFFT( unsigned flags, bool normalized )
486 {
487  doFFT( flags, FFTW_BACKWARD, normalized );
488 }
489 
490 /**
491  * Writes/Displays the object on an output stream.
492  * @param out the output stream where the object is written.
493  */
494 template <typename TSpace, typename T>
495 inline
496 void
497 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::selfDisplay ( std::ostream & out ) const
498 {
499  out << "[RealFFT]";
500 }
501 
502 /**
503  * Checks the validity/consistency of the object.
504  * @return 'true' if the object is valid, 'false' otherwise.
505  */
506 template <typename TSpace, typename T>
507 inline
508 bool
509 DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::isValid() const noexcept
510 {
511  return myStorage != nullptr;
512 }
513 
514 
515 
516 ///////////////////////////////////////////////////////////////////////////////
517 // Implementation of inline functions //
518 
519 template <
520  class TDomain,
521  typename T
522 >
523 inline
524 std::ostream&
525 DGtal::operator<< ( std::ostream & out,
526  const RealFFT<TDomain, T> & object )
527 {
528  object.selfDisplay( out );
529  return out;
530 }
531 
532 // //
533 ///////////////////////////////////////////////////////////////////////////////
534 
535