DGtal 1.4.0
Loading...
Searching...
No Matches
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.
50template <typename TSpace, typename T>
51inline
52DGtal::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.
71template <typename TSpace, typename T>
72inline
73DGtal::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
83template <typename TSpace, typename T>
84inline
85DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
86 ~RealFFT()
87{
88 FFTW::free( myStorage );
89}
90
91// Padding when using real datas.
92template <typename TSpace, typename T>
93inline
94std::size_t
95DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
96 getPadding() const noexcept
97{
98 return 2*myFreqExtent[0] - mySpatialExtent[0];
99}
100
101// Gets mutable spatial storage.
102template <typename TSpace, typename T>
103inline
104typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Real*
105DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
106 getSpatialStorage() noexcept
107{
108 return reinterpret_cast<Real*>(myStorage);
109}
110
111// Gets non-mutable spatial storage.
112template <typename TSpace, typename T>
113inline
114typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Real const*
115DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
116 getSpatialStorage() const noexcept
117{
118}
119
120// Gets mutable spatial image.
121template <typename TSpace, typename T>
122inline
123DGtal::ArrayImageAdapter<
124 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Real*,
125 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Domain
126>
127DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
128 getSpatialImage() noexcept
129{
130 return { getSpatialStorage(), myFullSpatialDomain, mySpatialDomain };
131}
132
133// Gets non-mutable spatial image.
134template <typename TSpace, typename T>
135inline
136DGtal::ArrayImageAdapter<
137 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Real const*,
138 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Domain
139>
140DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
141 getSpatialImage() const noexcept
142{
143 return { getSpatialStorage(), myFullSpatialDomain, mySpatialDomain };
144}
145
146// Gets mutable frequential raw storage.
147template <typename TSpace, typename T>
148inline
149typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Complex*
150DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
151 getFreqStorage() noexcept
152{
153 return reinterpret_cast<Complex*>(myStorage);
154}
155
156// Gets non-mutable frequential storage.
157template <typename TSpace, typename T>
158inline
159typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Complex const*
160DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
161 getFreqStorage() const noexcept
162{
163 return reinterpret_cast<Complex const*>(myStorage);
164}
165
166// Gets mutable frequential image.
167template <typename TSpace, typename T>
168inline
169DGtal::ArrayImageAdapter<
170 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Complex*,
171 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Domain
172>
173DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
174 getFreqImage() noexcept
175{
176 return { getFreqStorage(), getFreqDomain() };
177}
178
179// Get non-mutable frequential image.
180template <typename TSpace, typename T>
181inline
182DGtal::ArrayImageAdapter<
183 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Complex const*,
184 typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Domain
185>
186DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
187 getFreqImage() const noexcept
188{
189 return { getFreqStorage(), getFreqDomain() };
190}
191
192// Get spatial domain.
193template <typename TSpace, typename T>
194inline
195typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Domain const&
196DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
197 getSpatialDomain() const noexcept
198{
199 return mySpatialDomain;
200}
201
202// Get frequential domain.
203template <typename TSpace, typename T>
204inline
205typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Domain const&
206DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
207 getFreqDomain() const noexcept
208{
209 return myFreqDomain;
210}
211
212// Get spatial domain extent.
213template <typename TSpace, typename T>
214inline
215typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Point const&
216DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
217 getSpatialExtent() const noexcept
218{
219 return mySpatialExtent;
220}
221
222// Get frequential domain extent.
223template <typename TSpace, typename T>
224inline
225typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Point const&
226DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
227 getFreqExtent() const noexcept
228{
229 return myFreqExtent;
230}
231
232// Gets the extent of the scaled spatial domain.
233template <typename TSpace, typename T>
234inline
235typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::RealPoint
236DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::
237 getScaledSpatialExtent() const noexcept
238{
239 return myScaledSpatialExtent;
240}
241
242// Sets the extent of the scaled spatial domain.
243template <typename TSpace, typename T>
244inline
245void
246DGtal::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.
257template <typename TSpace, typename T>
258inline
259typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::RealPoint
260DGtal::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.
267template <typename TSpace, typename T>
268inline
269void
270DGtal::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.
277template <typename TSpace, typename T>
278inline
279typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::RealPoint
280DGtal::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.
293template <typename TSpace, typename T>
294inline
295typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Point
296DGtal::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.
311template <typename TSpace, typename T>
312inline
313typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::RealPoint
314DGtal::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.
331template <typename TSpace, typename T>
332inline
333typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Point
334DGtal::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.
363template <typename TSpace, typename T>
364inline
365typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Complex
366DGtal::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.
381template <typename TSpace, typename T>
382inline
383typename DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::Complex
384DGtal::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.
397template <typename TSpace, typename T>
398inline
399void
400DGtal::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.
416template <typename TSpace, typename T>
417inline
418void
419DGtal::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)
471template <typename TSpace, typename T>
472inline
473void
474DGtal::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)
481template <typename TSpace, typename T>
482inline
483void
484DGtal::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 */
494template <typename TSpace, typename T>
495inline
496void
497DGtal::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 */
506template <typename TSpace, typename T>
507inline
508bool
509DGtal::RealFFT<DGtal::HyperRectDomain<TSpace>, T>::isValid() const noexcept
510{
511 return myStorage != nullptr;
512}
513
514
515
516///////////////////////////////////////////////////////////////////////////////
517// Implementation of inline functions //
518
519template <
520 class TDomain,
521 typename T
522>
523inline
524std::ostream&
525DGtal::operator<< ( std::ostream & out,
526 const RealFFT<TDomain, T> & object )
527{
528 object.selfDisplay( out );
529 return out;
530}
531
532// //
533///////////////////////////////////////////////////////////////////////////////
534
535