SMIL  0.9.1
DFFT.hpp
1 /*
2  * Copyright (c) 2011-2016, Matthieu FAESSEL and ARMINES
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  *
8  * * Redistributions of source code must retain the above copyright
9  * notice, this list of conditions and the following disclaimer.
10  * * Redistributions in binary form must reproduce the above copyright
11  * notice, this list of conditions and the following disclaimer in the
12  * documentation and/or other materials provided with the distribution.
13  * * Neither the name of Matthieu FAESSEL, or ARMINES nor the
14  * names of its contributors may be used to endorse or promote products
15  * derived from this software without specific prior written permission.
16  *
17  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
18  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS AND CONTRIBUTORS BE
21  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27  * POSSIBILITY OF SUCH DAMAGE.
28  */
29 
30 #ifndef _D_FFT_HPP
31 #define _D_FFT_HPP
32 
33 #include "Core/include/private/DImage.hxx"
34 
35 #include <fftw3.h>
36 
37 namespace smil
38 {
60  template <class T1, class T2>
61  RES_T correlation(Image<T1> &imIn1, Image<T1> &imIn2, Image<T2> &imOut)
62  {
63  ASSERT_SAME_SIZE(&imIn1, &imIn2);
64 
65  size_t ncols = imIn1.getWidth();
66  // complex data column nbr
67  size_t nccols = ncols;
68  size_t nrows = imIn2.getHeight();
69 
70  imOut.setSize(ncols, nrows);
71 
72  size_t pixNbr = ncols * nrows;
73 
74  // Allocate arrays for FFT of src and tpl
75  double *src1_real = fftw_alloc_real(pixNbr);
76  fftw_complex *src1_compl = fftw_alloc_complex(nccols * nrows);
77  double *src2_real = fftw_alloc_real(pixNbr);
78  fftw_complex *src2_compl = fftw_alloc_complex(nccols * nrows);
79 
80  double *res_real = fftw_alloc_real(pixNbr);
81  fftw_complex *res_compl = fftw_alloc_complex(nccols * nrows);
82 
83  T1 *src1_data = imIn1.getPixels();
84  T1 *src2_data = imIn2.getPixels();
85  T2 *out_data = imOut.getPixels();
86 
87  // Copy image pixel values
88  for (size_t i = 0; i < pixNbr; i++) {
89  src1_real[i] = (double) src1_data[i];
90  src2_real[i] = (double) src2_data[i];
91  }
92 
93  // Create FFTW plans
94  fftw_plan forward1 = fftw_plan_dft_r2c_2d(nrows, ncols, src1_real,
95  src1_compl, FFTW_ESTIMATE);
96  fftw_plan forward2 = fftw_plan_dft_r2c_2d(nrows, ncols, src2_real,
97  src2_compl, FFTW_ESTIMATE);
98  fftw_plan backward = fftw_plan_dft_c2r_2d(nrows, ncols, res_compl, res_real,
99  FFTW_BACKWARD | FFTW_ESTIMATE);
100 
101  // Compute the FFT of the images
102  fftw_execute(forward1);
103  fftw_execute(forward2);
104 
105  // Compute the cross-correlation
106  for (size_t i = 0; i < pixNbr; i++) {
107  res_compl[i][0] = (src1_compl[i][0] * src2_compl[i][0] +
108  src1_compl[i][1] * src2_compl[i][1]);
109  res_compl[i][1] = (src1_compl[i][1] * src2_compl[i][0] -
110  src1_compl[i][0] * src2_compl[i][1]);
111 
112  double norm = sqrt(pow(res_compl[i][0], 2) + pow(res_compl[i][1], 2));
113  res_compl[i][0] /= norm;
114  res_compl[i][1] /= norm;
115  }
116 
117  // Compute result inverse fft
118  fftw_execute(backward);
119 
120  // Copy results to imOut
121  // Base results are between -1 and 1. We stretch/translate values to the
122  // output type value range.
123  for (size_t i = 0; i < pixNbr; i++) {
124  out_data[i] = ImDtTypes<T2>::min() +
125  (res_real[i] / pixNbr + 1) * ImDtTypes<T2>::cardinal() / 2;
126  }
127 
128  // Clear memory
129  fftw_destroy_plan(forward1);
130  fftw_destroy_plan(forward2);
131  fftw_destroy_plan(backward);
132  fftw_free(src1_real);
133  fftw_free(src1_compl);
134  fftw_free(src2_real);
135  fftw_free(src2_compl);
136  fftw_free(res_real);
137  fftw_free(res_compl);
138 
139  imOut.modified();
140 
141  return RES_OK;
142  }
143 
146 } // namespace smil
147 
148 #endif // _D_FFT_HPP
Author Vincent Morard Date : 7 march 2011 See : Morard V. and Dokladal P. and Decenciere E...
Definition: DApplyThreshold.hpp:36
RES_T correlation(Image< T1 > &imIn1, Image< T1 > &imIn2, Image< T2 > &imOut)
2D image (normalized) cross correlation using FFT.
Definition: DFFT.hpp:61
virtual RES_T setSize(size_t w, size_t h, size_t d=1, bool doAllocate=true)
Set the size of image.
Definition: DImage.hxx:319
size_t getHeight() const
Get image height.
Definition: DBaseImage.h:91
virtual void modified()
Trigger modified event (allows to force display update)
Definition: DImage.hxx:223
size_t getWidth() const
Get image width.
Definition: DBaseImage.h:87
Definition: DTypes.hpp:80
lineType getPixels() const
Get the pixels as a 1D array.
Definition: DImage.hpp:110