fft.hh 19 KB
Newer Older
Alexandre Abraham's avatar
Alexandre Abraham committed
1
// Copyright (C) 2007, 2008, 2009 EPITA Research and Development Laboratory
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
//
// This file is part of the Olena Library.  This library is free
// software; you can redistribute it and/or modify it under the terms
// of the GNU General Public License version 2 as published by the
// Free Software Foundation.
//
// This library 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 library; see the file COPYING.  If not, write to
// the Free Software Foundation, 51 Franklin Street, Fifth Floor,
// Boston, MA 02111-1307, USA.
//
// As a special exception, you may use this file as part of a free
// software library without restriction.  Specifically, if other files
// instantiate templates or use macros or inline functions from this
// file, or you compile this file and link it with other files to
// produce an executable, this file does not by itself cause the
// resulting executable to be covered by the GNU General Public
// License.  This exception does not however invalidate any other
// reasons why the executable file might be covered by the GNU General
// Public License.

#ifndef MLN_TRANSFORM_FFT_HH
# define MLN_TRANSFORM_FFT_HH

Alexandre Abraham's avatar
Alexandre Abraham committed
31
32
#  include <mln/core/image/image2d.hh>
#  include <mln/estim/min_max.hh>
33
#  include <mln/opt/at.hh>
34

Alexandre Abraham's avatar
Alexandre Abraham committed
35
36
37
#  include <complex>

#  include <fftw3.h>
38
39
40
41
42
43
44
45

namespace mln {

  namespace internal {

    /// Dispatch traits for fftw
    enum fft_dispatch { fft_cplx, fft_real };

Alexandre Abraham's avatar
Alexandre Abraham committed
46
47
48
    template <typename T>
    struct fft_trait;

49
50
51
    /*!
    ** \brief fft trait.
    */
Alexandre Abraham's avatar
Alexandre Abraham committed
52
53
    template <>
    struct fft_trait<double>
54
55
    {
      static const fft_dispatch		which = fft_real; ///< Real dispatch.
Alexandre Abraham's avatar
Alexandre Abraham committed
56
      typedef double			fftw_input; ///< Type of input.
57
58
59
60
61
62
63
    };

    /*!
    ** \brief fft trait.
    **
    ** \specialization for ntg::cplx<R, T>
    */
Alexandre Abraham's avatar
Alexandre Abraham committed
64
65
    template <typename T>
    struct fft_trait< std::complex<T> >
66
67
    {
      static const fft_dispatch		which = fft_cplx; ///< Complex dispatch.
Alexandre Abraham's avatar
Alexandre Abraham committed
68
      typedef std::complex <T>		fftw_input; ///< Type of input.
69
70
71
72
73
74
75
    };

    /*!
    ** _fft<ntg::cplx_representation, T>
    **
    ** \param T Data type.
    */
Alexandre Abraham's avatar
Alexandre Abraham committed
76
    template <class T>
77
78
79
80
81
82
83
84
    class _fft
    {
    public:
      /*!
      ** \brief Accessor to transformed image.
      **
      ** Const version.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
85
      const image2d< std::complex<T> > transformed_image() const
86
87
88
89
90
91
92
93
94
      {
	return trans_im;
      }

      /*!
      ** \brief Accessor to transformed image.
      **
      ** Non const version.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
95
      image2d< std::complex<T> >& transformed_image()
96
97
98
99
100
101
102
103
104
105
      {
	return trans_im;
      }

      /*!
      ** \brief Accessor to transformed image.
      **
      ** For each point p of the transformed image T, you get
      ** |T[p]|.
      **
Alexandre Abraham's avatar
Alexandre Abraham committed
106
      ** \param R Data type of the resulting image.
107
108
109
      **
      ** \arg ordered Kind of traversal.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
110
111
      template <class R>
      image2d<R> transformed_image_magn(bool ordered = true) const
112
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
113
	// FIXME : check that R is real
114

Alexandre Abraham's avatar
Alexandre Abraham committed
115
	image2d<R> new_im(trans_im.domain());
116
117

	if (ordered)
Alexandre Abraham's avatar
Alexandre Abraham committed
118
119
	  for (unsigned row = 0; row < new_im.nrows(); ++row)
	    for (unsigned col = 0; col < new_im.ncols(); ++col)
120
121
	      opt::at(new_im, row, col) =
		std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(),
Alexandre Abraham's avatar
Alexandre Abraham committed
122
				      (col + trans_im.ncols() / 2) % trans_im.ncols()));
123
124
	else
	  {
Alexandre Abraham's avatar
Alexandre Abraham committed
125
	    mln_piter(image2d< std::complex<T> >) it(trans_im.domain());
126
127

	    for_all(it)
Alexandre Abraham's avatar
Alexandre Abraham committed
128
	      new_im(it) = std::norm(trans_im(it));
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
	  }

	return new_im;
      }

      /*!
      ** \brief Accessor to transformed image.
      **
      ** For each point p of the transformed image T, you get
      ** |T[p]|.
      **
      ** \arg ordered Kind of traversal.
      */
      image2d<T> transformed_image_magn(bool ordered = true) const
      {
	return transformed_image_magn<T>(ordered);
      }

      /*!
      ** \brief Accessor to transformed image.
      **
      ** For each point p of the transformed image T, you get
      **  a clipped value of |T[p]|.\n
      **
Alexandre Abraham's avatar
Alexandre Abraham committed
153
      ** \param R Data type of the resulting image.
154
155
156
157
      **
      ** \arg clip Value used for clipping.
      ** \arg ordered Kind of traversal.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
158
159
160
161

      template <class R>
      image2d<R> transformed_image_clipped_magn(const double clip,
						bool ordered = true) const
162
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
163
	// Check that R is real
164

Alexandre Abraham's avatar
Alexandre Abraham committed
165
166
167
168
	image2d<R> new_im(trans_im.domain());
	// check that clip is >=0 and <=1 ?
	double max;
	mln_piter(image2d<T>) it(trans_im.domain());
169
170

	for_all(it)
Alexandre Abraham's avatar
Alexandre Abraham committed
171
172
	  if (std::norm(trans_im(it)) > max)
	    max = std::norm(trans_im(it));
173
174

	if (ordered)
Alexandre Abraham's avatar
Alexandre Abraham committed
175
176
	  for (unsigned row = 0; row < new_im.nrows(); ++row)
	    for (unsigned col = 0; col < new_im.ncols(); ++col)
177
	      {
178
		if (std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(),
Alexandre Abraham's avatar
Alexandre Abraham committed
179
180
					  (col + trans_im.ncols() / 2) %
					  trans_im.ncols())) >= max * clip)
181
		  opt::at(new_im, row, col) = mln_max(R);
182
		else
183
		  opt::at(new_im, row, col) =
Alexandre Abraham's avatar
Alexandre Abraham committed
184
		    (double) mln_max(R) *
185
		    std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(),
Alexandre Abraham's avatar
Alexandre Abraham committed
186
187
188
					  (col + trans_im.ncols() / 2) %
					  trans_im.ncols())) / (max * clip);
      }
189
	else
Alexandre Abraham's avatar
Alexandre Abraham committed
190
191
192
193
194
195
196
197
198
199
	  {
	    for_all(it)
	    {
	      if (std::norm(trans_im(it)) >= max * clip)
		new_im(it) = mln_max(R);
	      else
		new_im(it) = (double) mln_max(R) *
		  std::norm(trans_im(it)) / (max * clip);
	    }
	  }
200
201
202
203
204
205
206
207
208
209
210
211
212

	return new_im;
      }

      /*!
      ** \brief Accessor to transformed image.
      **
      ** For each point p of the transformed image T, you get
      **  a clipped value of |T[p]|.\n
      **
      ** \arg clip Value used for clipping.
      ** \arg ordered Kind of traversal.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
213
      image2d<T> transformed_image_clipped_magn(const double clip,
214
215
216
217
218
219
220
221
222
223
224
						bool ordered = true) const
      {
	return transformed_image_clipped_magn<T>(clip, ordered);
      }

      /*!
      ** \brief Accessor to transformed image.
      **
      ** For each point p of the transformed image T, you get
      ** a clipped value of |T[p]|.\n
      **
Alexandre Abraham's avatar
Alexandre Abraham committed
225
      ** \param R Data type of the resulting image.
226
227
228
      **
      ** \arg ordered Kind of traversal.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
229
230

      image2d<T> transformed_image_clipped_magn(bool ordered = true) const
231
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
232
	return transformed_image_clipped_magn<T>(1, ordered);
233
234
235
236
237
238
239
240
      }

      /*!
      ** \brief Accessor to transformed image.
      **
      ** For each point p of the transformed image T, you get
      ** a clipped value of |T[p]|.\n
      **
Alexandre Abraham's avatar
Alexandre Abraham committed
241
      ** \param R Data type of the resulting image.
242
243
244
      **
      ** \arg ordered Kind of traversal.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
245
246
247

      template <class R>
      image2d<R> transformed_image_clipped_magn(bool ordered = true) const
248
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
249
	return transformed_image_clipped_magn<R>(1, ordered);
250
251
252
253
254
255
256
257
258
259
260
261
262
      }

      // FIXME: Find a more elegant way to fix range interval on a and b.
      /*!
      ** \brief Accessor to transformed image.
      **
      ** For each point p of the transformed image T, you get
      ** a log translated value of |T[p]| on interval [a; b].\n
      **
      ** \arg a Lower bound.
      ** \arg b Upper bound.
      ** \arg ordered Kind of traversal.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
263
264
265
266
267

      template <class R>
      image2d<R> transformed_image_log_magn(double a,
					    double b,
					    bool ordered = true) const
268
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
269
270
271
272
273
	// Check that R is real
	// 0 <= a <= 1000
	// 0 <= b <= 1000

	image2d<R> new_im(trans_im.domain());
274

Alexandre Abraham's avatar
Alexandre Abraham committed
275
276
	double max = 0;
	mln_piter(image2d<R>) it(trans_im.domain());
277
278

	for_all(it)
Alexandre Abraham's avatar
Alexandre Abraham committed
279
280
281
	  if (std::norm(trans_im(it)) > max)
	    max = std::norm(trans_im(it));

282
283

	if (ordered)
Alexandre Abraham's avatar
Alexandre Abraham committed
284
285
	  for (unsigned row = 0; row < new_im.nrows(); ++row)
	    for (unsigned col = 0; col < new_im.ncols(); ++col)
286
287
	      opt::at(new_im, row, col) =
		log(a + b * std::norm(opt::at(trans_im, (row + trans_im.nrows() / 2) % trans_im.nrows(),
Alexandre Abraham's avatar
Alexandre Abraham committed
288
289
						  (col + trans_im.ncols() / 2) % trans_im.ncols()))) /
		log (a + b * max) * mln_max(R);
290
	else
Alexandre Abraham's avatar
Alexandre Abraham committed
291
292
293
294
295
296
297
	  {
	    mln_piter(image2d< std::complex<T> >) it(trans_im.domain());

	    for_all(it)
	      new_im(it) = log(a + b * std::norm(trans_im(it))) /
	      log (a + b * max) * mln_max(R);
	  }
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312

	return new_im;
      }

      // FIXME: Find a more elegant way to fix boundaries of a and b.
      /*!
      ** \brief Accessor to transformed image.
      **
      ** For each point p of the transformed image T, you get
      ** a log translated value of |T[p]| on interval [a; b].\n
      **
      ** \arg a Lower bound.
      ** \arg b Upper bound.
      ** \arg ordered Kind of traversal.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
313
314
315
316

	image2d<T> transformed_image_log_magn(double a,
					      double b,
					      bool ordered = true) const
317
318
319
320
321
322
323
324
325
326
      {
	return transformed_image_log_magn<T>(a, b, ordered);
      }

      /*!
      ** \brief Accessor to transformed image.
      **
      ** For each point p of the transformed image T, you get
      ** a log translated value of |T[p]| on interval [1; 100].\n
      **
Alexandre Abraham's avatar
Alexandre Abraham committed
327
      ** \param R Data type of the resulting image.
328
329
330
      **
      ** \arg ordered Kind of traversal.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
331
332
333

      template <class R>
      image2d<R> transformed_image_log_magn(bool ordered = true) const
334
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
335
	return transformed_image_log_magn<R>(1, 100, ordered);
336
337
338
339
340
341
342
343
344
345
      }

      /*!
      ** \brief Accessor to transformed image.
      **
      ** For each point p of the transformed image T, you get
      ** a log translated value of |T[p]| on interval [1; 100].\n
      **
      ** \arg ordered Kind of traversal.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
346

347
348
349
350
351
      image2d<T> transformed_image_log_magn(bool ordered = true) const
      {
	return transformed_image_log_magn<T>(1, 100, ordered);
      }

Alexandre Abraham's avatar
Alexandre Abraham committed
352

353
354
355
356
357
358
359
      /*!
      ** \brief Destructor.
      **
      ** Let memory free for other big images !!!
      */
      ~_fft()
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
360
361
362
363
	fftw_free(in);
	fftw_free(out);
	fftw_destroy_plan(p);
	fftw_destroy_plan(p_inv);
364
365
366
367
368
      }

    protected:

      typename fft_trait<T>::fftw_input	*in; ///< Input image.
Alexandre Abraham's avatar
Alexandre Abraham committed
369
370
371
372
      std::complex<T>			*out; ///< Complex image.
      fftw_plan				p; ///< Plan.
      fftw_plan				p_inv; ///< inverted plan.
      image2d< std::complex<T> >	trans_im; ///< Transformed image.
373
374
375
376
377

    };

  } // end of namespace internal

Alexandre Abraham's avatar
Alexandre Abraham committed
378
  namespace transform {
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394

    /*!
    ** \brief fft class declaration.
    **
    ** \param T Data type.
    ** \param which Dispatch.
    */
    template <class T,
	      internal::fft_dispatch which = internal::fft_trait<T>::which >
    class fft;

    /*!
    ** \brief fft specialization for fft real dispatch.
    **
    ** \param T Data type.
    */
Alexandre Abraham's avatar
Alexandre Abraham committed
395
396
    template <class T>
    class fft<T, internal::fft_real> : public internal::_fft<T>
397
398
399
400
401
402
403
404
405
406
407
    {

    public:

      /*!
      ** \brief Constructor.
      **
      ** Initialization of data in order to compute the fft.
      **
      ** \arg original_im Image to process.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
408
409
      template <typename D>
      fft(const image2d<D>& original_im)
410
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
411
412
413
	this->in  = (T*) fftw_malloc(original_im.nrows() * original_im.ncols() * sizeof(T));
	this->out = (std::complex<T>*)
	  fftw_malloc(original_im.nrows() * (original_im.ncols() / 2 + 1) * sizeof(std::complex<T>));
414

Alexandre Abraham's avatar
Alexandre Abraham committed
415
416
	for (unsigned row = 0; row < original_im.nrows(); ++row)
	  for (unsigned col = 0; col < original_im.ncols(); ++col)
417
	    this->in[row * original_im.ncols() + col] = opt::at(original_im, row, col);
418

Alexandre Abraham's avatar
Alexandre Abraham committed
419
420
421
422
	this->p = fftw_plan_dft_r2c_2d (original_im.nrows(), original_im.ncols(),
					this->in, reinterpret_cast<fftw_complex*>(this->out), FFTW_ESTIMATE);
	this->p_inv = fftw_plan_dft_r2c_2d (original_im.nrows(), original_im.ncols(),
					    this->in, reinterpret_cast<fftw_complex*>(this->out), FFTW_ESTIMATE);
423

Alexandre Abraham's avatar
Alexandre Abraham committed
424
	this->trans_im = image2d< std::complex<T> >(original_im.domain());
425
426
427
428
429
      }

      /*!
      ** \brief Compute and return the transform.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
430
      image2d< std::complex<T> > transform()
431
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
432
	fftw_execute(this->p);
433
434
435

	unsigned denom = this->trans_im.nrows() * this->trans_im.ncols();
	int i = 0;
Alexandre Abraham's avatar
Alexandre Abraham committed
436
437
	for (unsigned row = 0; row < this->trans_im.nrows(); ++row)
	  for (unsigned col = 0; col <= this->trans_im.ncols() / 2; ++col)
438
	    {
Alexandre Abraham's avatar
Alexandre Abraham committed
439
	      this->out[i] = std::complex<T> (this->out[i].real() / denom, this->out[i].imag() / denom);
440
	      opt::at(this->trans_im, row, col) = this->out[i];
441
442
	      ++i;
	    }
Alexandre Abraham's avatar
Alexandre Abraham committed
443
444
	for (unsigned row = 0; row < this->trans_im.nrows(); ++row)
	  for (unsigned col = this->trans_im.ncols() - 1; col > this->trans_im.ncols() / 2; --col)
Alexandre Abraham's avatar
Alexandre Abraham committed
445
	    opt::at(this->trans_im, row, col) = opt::at(this->trans_im, this->trans_im.nrows() - row - 1,
Alexandre Abraham's avatar
Alexandre Abraham committed
446
							    this->trans_im.ncols() - col - 1);
447
448
449
450
451
452
	return this->trans_im;
      }

      /*!
      ** \brief Compute and return the invert transform.
      **
Alexandre Abraham's avatar
Alexandre Abraham committed
453
      ** \param R Data type of output image.
454
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
455
456
      template <class R>
      image2d<R> transform_inv()
457
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
458
	// FIXME : Check that R is real
459

Alexandre Abraham's avatar
Alexandre Abraham committed
460
461
462
	for (unsigned row = 0; row < this->trans_im.nrows(); ++row)
	  for (unsigned col = 0; col <= this->trans_im.ncols() / 2; ++col)
	    this->out[row * (this->trans_im.ncols() / 2 + 1) + col] =
463
	      opt::at(this->trans_im, row, col).real();
Alexandre Abraham's avatar
Alexandre Abraham committed
464
465

	fftw_execute(this->p_inv);
466

Alexandre Abraham's avatar
Alexandre Abraham committed
467
	image2d<R> new_im(this->trans_im.domain());
468
	int i = 0;
Alexandre Abraham's avatar
Alexandre Abraham committed
469
470
	for (unsigned row = 0; row < this->trans_im.nrows(); ++row)
	  for (unsigned col = 0; col < this->trans_im.ncols(); ++col)
471
	    {
472
	      opt::at(new_im, row, col) = (this->in[i] >= mln_min(R) ?
Alexandre Abraham's avatar
Alexandre Abraham committed
473
474
475
476
				     (this->in[i] <= mln_max(R) ?
				      (R)this->in [i] :
				      mln_min(R)) :
				     mln_max(R));
477
478
479
480
481
482
483
484
485
	      ++i;
	    }
	return new_im;
      }

      /*!
      ** \brief Shift zero-frequency component of discrete Fourier transform
      ** to center of spectrum.
      **
Alexandre Abraham's avatar
Alexandre Abraham committed
486
      ** \param R The data type of the image returned.
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
      **
      ** The zero-frequency component of discrete Fourier transform are moved
      ** to center of the image :
      **
      ** \htmlonly
      ** <table>
      **   <tr><td>1</td><td>2</td></tr>
      **   <tr><td>3</td><td>4</td></tr>
      ** </table>
      ** becomes
      ** <table>
      **   <tr><td>4</td><td>3</td></tr>
      **   <tr><td>2</td><td>1</td></tr>
      ** </table>
      ** \endhtmlonly
      **
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
504
505
506
507

      /*
      template <class R>
      image2d<R> shift_transform_inv()
508
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
509
510
	image2d<R> t = transform_inv<R>();
	image2d<R> st(t.size());
511
512
513

	// We have to exchange t_1 with t_1_dest and not directly t_3 because
	// they have not he same size.
Alexandre Abraham's avatar
Alexandre Abraham committed
514
	typedef morpher::piece_morpher< image2d<R> > piece_t;
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
	piece_t t_1(t, dpoint2d(0, 0),
		    image2d_size((t.size().nrows() - 1) / 2,
				 (t.size().ncols() - 1) / 2,
				 t.border()));
	piece_t t_1_dest(st, dpoint2d(t.nrows() - t_1.nrows(),
				     t.ncols() - t_1.ncols()),
			 image2d_size(t_1.nrows(), t_1.ncols(),
				      t.border()));
	piece_t t_2(t, dpoint2d(0, (t.size().ncols() - 1) / 2),
		    image2d_size((t.size().nrows() - 1) / 2,
				 t.size().ncols() - (t.size().ncols() - 1) / 2,
				 t.border()));
	piece_t t_2_dest(st, dpoint2d(t.nrows() - t_2.nrows(), 0),
			 image2d_size(t_2.nrows(), t_2.ncols(),
				      t.border()));
	piece_t t_3(t, dpoint2d((t.size().nrows() - 1) / 2, 0),
		    image2d_size(t.size().nrows() - (t.size().nrows() - 1) / 2,
				 (t.size().ncols() - 1) / 2,
				 t.border()));
	piece_t t_3_dest(st, dpoint2d(0, t.ncols() - t_3.ncols()),
			 image2d_size(t_3.nrows(), t_3.ncols(),
				      t.border()));
	piece_t t_4(t, dpoint2d((t.size().nrows() - 1) / 2,
				(t.size().ncols() - 1) / 2),
		    image2d_size(t.size().nrows() - (t.size().nrows() - 1) / 2,
				 t.size().ncols() - (t.size().ncols() - 1) / 2,
				 t.border()));
	piece_t t_4_dest(st, dpoint2d(0, 0),
			 image2d_size(t_4.nrows(), t_4.ncols(),
				      t.border()));

	oln_iter_type(piece_t) i1(t_1);
	for_all(i1)
	  t_1_dest[i1] = t_1[i1];
	oln_iter_type(piece_t) i2(t_2);
	for_all(i2)
	  t_2_dest[i2] = t_2[i2];
	oln_iter_type(piece_t) i3(t_3);
	for_all(i3)
	  t_3_dest[i3] = t_3[i3];
	oln_iter_type(piece_t) i4(t_4);
	for_all(i4)
	  t_4_dest[i4] = t_4[i4];

	return st;
      }
Alexandre Abraham's avatar
Alexandre Abraham committed
561
562
      */

563
564
565
566
567
568
569
570
571
572
573
574
575

      /*!
      ** \brief Compute and return the invert transform.
      */
      image2d<T> transform_inv()
      {
	return transform_inv<T>();
      }

      /*!
      ** \brief Shift zero-frequency component of discrete Fourier transform
      ** to center of spectrum.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
576
      /*
577
578
579
580
      image2d<T> shift_transform_inv()
      {
	return shift_transform_inv<T>();
      }
Alexandre Abraham's avatar
Alexandre Abraham committed
581
      */
582
583
584
585
586
587
588
589
590

    };

    /*!
    ** \brief fft specialization for fft complex dispatch.
    **
    ** \param T Data type.
    ** \param R Complex representation kind.
    */
Alexandre Abraham's avatar
Alexandre Abraham committed
591
592
    template <class T>
    class fft<T, internal::fft_cplx> : public internal::_fft<T>
593
594
595
596
597
598
599
600
601
602
603
    {

    public:

      /*!
      ** \brief Constructor.
      **
      ** Initialization of data in order to compute the fft.
      **
      ** \arg original_im Image to process.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
604
      fft(const image2d< std::complex<T> >& original_im)
605
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
606
607
608
609
	this->in = fftw_malloc(sizeof(std::complex<T>) *
			       original_im.nrows() * original_im.ncols());
	this->out = fftw_malloc(sizeof(std::complex<T>) *
				original_im.nrows() * original_im.ncols());
610

Alexandre Abraham's avatar
Alexandre Abraham committed
611
612
	for (unsigned row = 0; row < original_im.nrows(); ++row)
	  for (unsigned col = 0; col < original_im.ncols(); ++col)
613
614
615
616
617
618
619
	    {
	      this->in[row * original_im.ncols() + col].re =
		original_im(row, col).real();
	      this->in[row * original_im.ncols() + col].im =
		original_im(row, col).imag();
	    }

Alexandre Abraham's avatar
Alexandre Abraham committed
620
621
622
623
624
	this->p = fftw_plan_dft_2d(original_im.nrows(), original_im.ncols(),
				   this->in, this->out,
				   FFTW_FORWARD, FFTW_ESTIMATE);
	this->p_inv = fftw2d_plan_dft_2d(original_im.nrows(), original_im.ncols(),
					 this->in, this->out,
625
626
					 FFTW_BACKWARD, FFTW_ESTIMATE);

Alexandre Abraham's avatar
Alexandre Abraham committed
627
	this->trans_im = image2d< std::complex<T> >(original_im.domain());
628
629
630
631
632
      }

      /*!
      ** \brief Compute and return the transform.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
633
      image2d< std::complex<T> > transform()
634
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
635
	fftw_execute(this->p);
636
637
638

	unsigned denom = this->trans_im.nrows() * this->trans_im.ncols();
	int i = 0;
Alexandre Abraham's avatar
Alexandre Abraham committed
639
640
	for (unsigned row = 0; row < this->trans_im.nrows(); ++row)
	  for (unsigned col = 0; col < this->trans_im.ncols(); ++col)
641
642
643
	    {
	      this->out[i].re /= denom;
	      this->out[i].im /= denom;
Alexandre Abraham's avatar
Alexandre Abraham committed
644
645
	      this->trans_im(row, col) = std::complex<double>(this->out[i].re,
							 this->out[i].im);
646
647
648
649
650
651
652
653
	      ++i;
	    }
	return this->trans_im;
      }

      /*!
      ** \brief Compute and return the invert transform.
      **
Alexandre Abraham's avatar
Alexandre Abraham committed
654
      ** \param R Data type of output image.
655
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
656
657
      template <class R>
      image2d< std::complex<R> > transform_inv()
658
      {
Alexandre Abraham's avatar
Alexandre Abraham committed
659
660
	for (unsigned row = 0; row < this->trans_im.nrows(); ++row)
	  for (unsigned col = 0; col < this->trans_im.ncols(); ++col)
661
662
663
664
665
	    {
	      this->out[row * this->trans_im.ncols() + col].re = this->trans_im(row, col).real();
	      this->out[row * this->trans_im.ncols() + col].im = this->trans_im(row, col).imag();
	    }

Alexandre Abraham's avatar
Alexandre Abraham committed
666
667
668
	fftw_execute(this->p_inv);

	image2d< std::complex<R> > new_im(this->trans_im.size());
669
	int i = 0;
Alexandre Abraham's avatar
Alexandre Abraham committed
670
671
	for (unsigned row = 0; row < this->trans_im.nrows(); ++row)
	  for (unsigned col = 0; col < this->trans_im.ncols(); ++col)
672
673
674
675
676
677
678
679
680
681
	    {
	      new_im(row, col) = this->in[i];
	      ++i;
	    }
	return new_im;
      }

      /*!
      ** \brief Compute and return the invert transform.
      */
Alexandre Abraham's avatar
Alexandre Abraham committed
682
      image2d< std::complex<T> > transform_inv()
683
684
685
686
687
688
      {
	return transform_inv<T>();
      }

    };

Alexandre Abraham's avatar
Alexandre Abraham committed
689
  } // end of namespace transform
690
691
692
693

} // end of namespace oln

#endif // ! MLN_TRANSFORM_FFT_HH