Line data Source code
1 : /** \file mftT.hpp
2 : * \brief The Matrix Fourier Transform interface
3 : * \ingroup ft_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2024-2025 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef mftT_hpp
28 : #define mftT_hpp
29 :
30 : #pragma GCC system_header
31 : #include <Eigen/Dense>
32 :
33 : #include "fftwTemplates.hpp"
34 : #include "ftTypes.hpp"
35 : #include "../../math/constants.hpp"
36 :
37 : /** \addtogroup mft
38 : *
39 : * The Matrix Fourier Transform is the <a href="https://en.wikipedia.org/wiki/DFT_matrix">
40 : * matrix-multiplication implementation of the
41 : * Discrete Fourier Transform</a>. It is slower than the FFT given the same size matrix, e.g.
42 : * <a href="https://pages.hmc.edu/ruye/e161/lectures/fourier/node11.html">in 2D</a>
43 : * \f$O(N^3)\f$ vs \f$O(N^2\ln(N))\f$. However, compared to zero-padding,
44 : * it can provide advantages in both speed
45 : * and memory required for problems where oversampling is needed but only over a small region
46 : * of the output.
47 : */
48 : namespace mx
49 : {
50 : namespace math
51 : {
52 : namespace ft
53 : {
54 :
55 : template <typename _inputT, typename _outputT, size_t _rank, int _cudaGPU = 0>
56 : class mftT;
57 :
58 : /// Matrix Fourier Transforms
59 : /** Calculates the Discrete Fourier Transform (DFT) using matrix multiplication. This
60 : * is normally much less efficient than the FFT, but for large oversampling (padding)
61 : * the Matrix FT (MFT) will be much more space efficient and faster at the cost of
62 : * "field of view".
63 : *
64 : * This interface is modeled after the \ref fftT<_inputT, _outputT, _rank, 0> "fftT" interface to fftw, and is interoperable with it.
65 : * That means that using, e.g., mftT (unshifted, not-oversampled) for the forward transform
66 : * and \ref fftT<_inputT, _outputT, _rank, 0> "fftT" for the backward
67 : * transform is equivalent to using one or the other for both. Note that
68 : * the MFT is normalized as in fftw, so by 1/N on the forward and by 1 on the backward.
69 : *
70 : * Oversampling is optionally included in the transform, which is equivalent to zero padding
71 : * in terms of increased resolution \cite soummer_2007.
72 : * The cost is that the output domain is truncated by the oversampling factor. I.e. if we
73 : * oversample by a factor of 10, only 1/10th of the output transform is available.
74 : *
75 : * The output can be shifted as part of the MFT calculation, which is implemented similarly
76 : * to <a href="https://www.mathworks.com/matlabcentral/fileexchange/18401-efficient-subpixel-image-registration-by-cross-correlation">
77 : * this matlab code</a>.
78 : *
79 : * Note that when either oversampling or shifting is done on a forward (backward) transform,
80 : * a subsequent backward (forward) transform will not in general be the inverse.
81 : *
82 : * \tparam inputT is the input type of the transform, only complex types are suppored by mftT
83 : * \tparam outputT is the output type of the transform, only complex types are suppored by mftT
84 : * \tparam _rank is the rank of the transform. Currently only rank 2 is implemented.
85 : *
86 : * \ingroup mft
87 : */
88 : template <typename _inputT, typename _outputT, size_t _rank>
89 : class mftT<_inputT, _outputT, _rank, 0>
90 : {
91 : typedef _inputT inputT;
92 : typedef _outputT outputT;
93 :
94 : typedef typename fftwTypeSpec<inputT, outputT>::realT realT;
95 :
96 : typedef typename fftwTypeSpec<_inputT, _outputT>::complexT complexT;
97 :
98 : static const size_t rank = _rank;
99 :
100 : typedef Eigen::Array<inputT, -1, -1> eigenArrayInT;
101 : typedef Eigen::Array<outputT, -1, -1> eigenArrayOutT;
102 :
103 : typedef Eigen::Matrix<complexT, -1, -1> eigenMatrixT;
104 :
105 : protected:
106 : dir m_dir{ dir::forward }; /**< Direction of this MFT, either dir::forward (default)
107 : or dir::backward */
108 :
109 : int m_szX{ 0 }; ///< Size of the x dimension
110 : int m_szY{ 0 }; ///< Size of the y dimension
111 : int m_szZ{ 0 }; ///< size of the z dimension
112 :
113 : float m_osFac{ 1 }; ///< The oversampling factor
114 :
115 : realT m_xOff{ 0 }; ///< The offset in the rows direction for the center of the DFT.
116 : realT m_yOff{ 0 }; ///< The offset in the columns direction for the center of the DFT.
117 : public:
118 : eigenMatrixT m_dftR; ///< DFT matrix for the rows
119 : eigenMatrixT m_dftC; ///< DFT matrix for the columnss
120 :
121 : public:
122 : /// Default c'tor
123 : mftT();
124 :
125 : /// Constructor for rank 1 MFT.
126 : template <size_t crank = _rank>
127 : mftT( int nx, ///< [in] the desired size of the MFT
128 : dir ndir = dir::forward, /**< [in] [optional] direction of this MFT, either dir::forward
129 : (default) or dir::backward */
130 : realT xOff = 0, /**< [in] [optional] the x offset of the center of the
131 : transformed array. Default 0.*/
132 : realT osFac = 1.0, /**< [in] [optional] the oversampling factor. Default 1. */
133 : typename std::enable_if<crank == 1>::type * = 0 ) = delete;
134 :
135 : /// Constructor for rank 2 MFT.
136 : template <size_t crank = _rank>
137 : mftT( int nx, ///< [in] the desired x size of the MFT
138 : int ny, ///< [in] the desired y size of the MFT
139 : dir ndir = dir::forward, /**< [in] [optional] direction of this MFT, either dir::forward
140 : (default) or dir::backward */
141 : realT xOff = 0, /**< [in] [optional] the x offset of the center of the
142 : transformed array. Default 0.*/
143 : realT yOff = 0, /**< [in] [optional] the y offset of the center of the
144 : transformed array Default 0.*/
145 : realT osFac = 1.0, /**< [in] [optional] the oversampling factor. Default 1. */
146 : typename std::enable_if<crank == 2>::type * = 0 );
147 :
148 : /// Constructor for rank 3 MFT.
149 : template <size_t crank = _rank>
150 : mftT( int nx, ///< [in] the desired x size of the MFT
151 : int ny, ///< [in] the desired y size of the MFT
152 : int nz, ///< [in] the desired z size of the MFT
153 : dir ndir = dir::forward, /**< [in] [optional] direction of this MFT, either dir::forward
154 : (default) or dir::backward */
155 : realT xOff = 0, /**< [in] [optional] the x offset of the center of the
156 : transformed array. Default 0.*/
157 : realT yOff = 0, /**< [in] [optional] the y offset of the center of the
158 : transformed array Default 0.*/
159 : realT zOff = 0, /**< [in] [optional] the z offset of the center of the
160 : transformed array Default 0.*/
161 : realT osFac = 1.0, /**< [in] [optional] the oversampling factor. Default 1. */
162 : typename std::enable_if<crank == 3>::type * = 0 ) = delete;
163 :
164 : /// Planning routine for rank 2 transforms.
165 : template <size_t crank = _rank>
166 : void plan( int nx, ///< [in] the desired x size of the MFT
167 : int ny, ///< [in] the desired y size of the MFT
168 : dir ndir = dir::forward, /**< [in] [optional] direction of this MFT, either dir::forward
169 : (default) or dir::backward */
170 : realT xOff = 0, /**< [in] [optional] the x offset of the center of the
171 : transformed array. Default 0.*/
172 : realT yOff = 0, /**< [in] [optional] the y offset of the center of the
173 : transformed array Default 0.*/
174 : realT osFac = 1.0, /**< [in] [optional] the oversampling factor. Default 1. */
175 : typename std::enable_if<crank == 2>::type * = 0 );
176 :
177 : /// Conduct the MFT
178 : void operator()( eigenArrayOutT &out, /**< [out] the output of the DFT */
179 : const eigenArrayInT &in /**< [in] the input to the DFT */ ) const;
180 : };
181 :
182 : template <typename inputT, typename outputT, size_t rank>
183 12 : mftT<inputT, outputT, rank, 0>::mftT()
184 : {
185 12 : }
186 :
187 : template <typename inputT, typename outputT, size_t rank>
188 : template <size_t crank>
189 : mftT<inputT, outputT, rank, 0>::mftT(
190 : int nx, int ny, dir ndir, realT xoff, realT yoff, realT osFac, typename std::enable_if<crank == 2>::type * )
191 : {
192 : plan( nx, ny, ndir, xoff, yoff, osFac );
193 : }
194 :
195 : template <typename inputT, typename outputT, size_t rank>
196 : template <size_t crank>
197 0 : void mftT<inputT, outputT, rank, 0>::plan(
198 : int nx, int ny, dir ndir, realT xOff, realT yOff, realT osFac, typename std::enable_if<crank == 2>::type * )
199 : {
200 0 : if(m_szX == nx && m_szY == ny && m_dir == ndir && m_xOff == xOff && m_yOff == yOff && m_osFac == osFac)
201 : {
202 0 : return;
203 : }
204 :
205 0 : m_szX = nx;
206 0 : m_szY = ny;
207 0 : m_dir = ndir;
208 0 : m_xOff = xOff;
209 0 : m_yOff = yOff;
210 0 : m_osFac = osFac;
211 :
212 0 : if( m_szX != m_szY )
213 : {
214 0 : throw std::invalid_argument( "MFT of non-square size is not implemented. nx must equal ny." );
215 : }
216 :
217 : // These should depend on szX too.
218 0 : m_dftR.resize( m_szX, m_szX );
219 0 : m_dftC.resize( m_szX, m_szX );
220 :
221 : // There is probably an osNx and an osNy?
222 0 : realT osN = m_szX * m_osFac;
223 :
224 0 : if( m_dir == dir::forward )
225 : {
226 0 : realT sign = -1;
227 0 : realT norm = 1.0 / ( m_szX * m_szY );
228 :
229 0 : for( int cc = 0; cc < m_szY; ++cc )
230 : {
231 0 : realT ccx = cc;
232 :
233 0 : for( int rr = 0; rr < m_szX; ++rr )
234 : {
235 0 : realT x = ( rr - m_xOff ) * ccx / osN;
236 :
237 0 : m_dftR( rr, cc ) = norm * exp( complexT( { 0, sign * 2 * pi<realT>() * x } ) );
238 :
239 0 : realT rrx = rr;
240 :
241 0 : x = rrx * ( cc - m_yOff ) / osN;
242 :
243 0 : m_dftC( rr, cc ) = norm * exp( complexT( { 0, sign * 2 * pi<realT>() * x } ) );
244 : }
245 : }
246 : }
247 : else
248 : {
249 0 : realT sign = +1;
250 0 : realT norm = 1.0;
251 :
252 0 : for( int cc = 0; cc < m_szY; ++cc )
253 : {
254 0 : realT ccx = cc;
255 0 : if( ccx > m_szY / 2 )
256 0 : ccx = -1 * ( m_szY - ccx );
257 :
258 0 : for( int rr = 0; rr < m_szX; ++rr )
259 : {
260 0 : realT rrx = rr;
261 0 : if( rrx > m_szX / 2 )
262 0 : rrx = -1 * ( m_szX - rrx );
263 :
264 0 : realT x = rrx * ( cc - m_xOff ) / osN;
265 :
266 0 : m_dftR( cc, rr ) = norm * exp( complexT( { 0, sign * 2 * pi<realT>() * x } ) );
267 :
268 0 : x = ( rr - m_yOff ) * ccx / osN;
269 :
270 0 : m_dftC( cc, rr ) = norm * exp( complexT( { 0, sign * 2 * pi<realT>() * x } ) );
271 : }
272 : }
273 : }
274 : }
275 :
276 : template <typename inputT, typename outputT, size_t rank>
277 0 : void mftT<inputT, outputT, rank, 0>::operator()( eigenArrayOutT &out, const eigenArrayInT &in ) const
278 : {
279 0 : out = ( m_dftR * in.matrix() * m_dftC ).array();
280 0 : }
281 :
282 : } // namespace ft
283 : } // namespace math
284 : } // namespace mx
285 :
286 : #endif // mdft_hpp
|
