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'' AND ANY
18  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
19  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
20  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS AND CONTRIBUTORS BE LIABLE FOR ANY
21  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
22  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
23  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
24  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27  */
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 
38 namespace smil
39 {
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  size_t nccols = ncols; // complex data column nbr
67  size_t nrows = imIn2.getHeight();
68 
69  imOut.setSize(ncols, nrows);
70 
71  size_t pixNbr = ncols * nrows;
72 
73  //Allocate arrays for FFT of src and tpl
74  double *src1_real = fftw_alloc_real(pixNbr);
75  fftw_complex *src1_compl = fftw_alloc_complex(nccols*nrows);
76  double *src2_real = fftw_alloc_real(pixNbr);
77  fftw_complex *src2_compl = fftw_alloc_complex(nccols*nrows);
78 
79  double *res_real = fftw_alloc_real(pixNbr);
80  fftw_complex *res_compl = fftw_alloc_complex(nccols*nrows);
81 
82  T1 *src1_data = imIn1.getPixels();
83  T1 *src2_data = imIn2.getPixels();
84  T2 *out_data = imOut.getPixels();
85 
86  // Copy image pixel values
87  for(size_t i = 0; i < pixNbr; i++)
88  {
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, src1_compl, FFTW_ESTIMATE);
95  fftw_plan forward2 = fftw_plan_dft_r2c_2d(nrows, ncols, src2_real, src2_compl, FFTW_ESTIMATE);
96  fftw_plan backward = fftw_plan_dft_c2r_2d(nrows, ncols, res_compl, res_real, FFTW_BACKWARD | FFTW_ESTIMATE);
97 
98  // Compute the FFT of the images
99  fftw_execute(forward1);
100  fftw_execute(forward2);
101 
102  // Compute the cross-correlation
103  for(size_t i = 0; i < pixNbr ; i++ )
104  {
105  res_compl[i][0] = (src1_compl[i][0]*src2_compl[i][0] + src1_compl[i][1]*src2_compl[i][1]);
106  res_compl[i][1] = (src1_compl[i][1]*src2_compl[i][0] - src1_compl[i][0]*src2_compl[i][1]);
107 
108  double norm = sqrt( pow(res_compl[i][0], 2) + pow(res_compl[i][1], 2) );
109  res_compl[i][0] /= norm;
110  res_compl[i][1] /= norm;
111  }
112 
113  // Compute result inverse fft
114  fftw_execute(backward);
115 
116  // Copy results to imOut
117  // Base results are between -1 and 1. We stretch/translate values to the output type value range.
118  for(size_t i=0;i<pixNbr;i++)
119  {
120  out_data[i] = ImDtTypes<T2>::min() + ( res_real[i]/pixNbr + 1) * ImDtTypes<T2>::cardinal()/2;
121  }
122 
123  // Clear memory
124  fftw_destroy_plan(forward1);
125  fftw_destroy_plan(forward2);
126  fftw_destroy_plan(backward);
127  fftw_free(src1_real);
128  fftw_free(src1_compl);
129  fftw_free(src2_real);
130  fftw_free(src2_compl);
131  fftw_free(res_real);
132  fftw_free(res_compl);
133 
134  imOut.modified();
135 
136  return RES_OK;
137  }
138 
141 } // namespace smil
142 
143 #endif // _D_FFT_HPP
Definition: DColorConvert.h:38
lineType getPixels() const
Get the pixels as a 1D array.
Definition: DImage.hpp:110
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
RES_T correlation(Image< T1 > &imIn1, Image< T1 > &imIn2, Image< T2 > &imOut)
2D image (normalized) cross correlation using FFT.
Definition: DFFT.hpp:61
Definition: DTypes.hpp:78