diff --git a/Base/Math/FourierTransform.cpp b/Base/Math/FourierTransform.cpp index f69711a59f7b44113e6fba26c7a5d5cbe24bb906..b5d460c5e8816482cbdca9d9c05537461f65ed19 100644 --- a/Base/Math/FourierTransform.cpp +++ b/Base/Math/FourierTransform.cpp @@ -22,48 +22,30 @@ namespace { template <typename T> void fftshift_1d(std::vector<T>& result, bool inverse) { // Centering FT around zero frequency - size_t col_shift; - if (result.size() % 2 == 0) - col_shift = result.size() / 2; - else - col_shift = result.size() / 2 + 1; + int col_shift = (result.size() + 1) / 2; if (inverse) - std::rotate(result.rbegin(), result.rbegin() + static_cast<int>(col_shift), result.rend()); + std::rotate(result.rbegin(), result.rbegin() + col_shift, result.rend()); else - std::rotate(result.begin(), result.begin() + static_cast<int>(col_shift), result.end()); + std::rotate(result.begin(), result.begin() + col_shift, result.end()); } template <typename T> void fftshift_2d(std::vector<std::vector<T>>& result, bool inverse) { // Centering FT around zero frequency - size_t row_shift; - if (result.size() % 2 == 0) - row_shift = result.size() / 2; - else - row_shift = result.size() / 2 + 1; - - size_t col_shift; - if (result[0].size() % 2 == 0) - col_shift = result[0].size() / 2; - else - col_shift = result[0].size() / 2 + 1; - - // First, shifting the rows - if (inverse) - std::rotate(result.rbegin(), result.rbegin() + static_cast<int>(row_shift), result.rend()); - else - std::rotate(result.begin(), result.begin() + static_cast<int>(row_shift), result.end()); - - // Second, shifting the cols - if (inverse) - for (size_t i = 0; i < result.size(); i++) - std::rotate(result[i].rbegin(), result[i].rbegin() + static_cast<int>(col_shift), - result[i].rend()); - else - for (size_t i = 0; i < result.size(); i++) - std::rotate(result[i].begin(), result[i].begin() + static_cast<int>(col_shift), - result[i].end()); + int row_shift = (result.size() + 1) / 2; + int col_shift = (result[0].size() + 1) / 2; + + // First shift the rows, then the cols + if (inverse) { + std::rotate(result.rbegin(), result.rbegin() + row_shift, result.rend()); + for (auto& row : result) + std::rotate(row.rbegin(), row.rbegin() + col_shift, row.rend()); + } else { + std::rotate(result.begin(), result.begin() + row_shift, result.end()); + for (auto& row : result) + std::rotate(row.begin(), row.begin() + col_shift, row.end()); + } } } // namespace @@ -111,8 +93,8 @@ void FourierTransform::Workspace::clear() void FourierTransform::rfft(const double2d_t& source, complex2d_t& result) { // rows (h) and columns (w) of the input 'source' - int h = static_cast<int>(source.size()); - int w_real = static_cast<int>((!source.empty() ? source[0].size() : 0)); + int h = source.size(); + int w_real = !source.empty() ? source[0].size() : 0; // initialization init_r2c(h, w_real); @@ -127,7 +109,7 @@ void FourierTransform::rfft(const double2d_t& source, complex2d_t& result) void FourierTransform::irfft(const complex2d_t& source, double2d_t& result, int w_real) { // rows (h) of the input 'source' - int h = static_cast<int>(source.size()); + int h = source.size(); // initialization init_c2r(h, w_real); @@ -246,8 +228,8 @@ void FourierTransform::init(int h, int w_real) ws.w_real = w_real; ws.w_fftw = w_real / 2 + 1; - ws.arr_real = fftw_alloc_real(static_cast<size_t>(ws.h * ws.w_real)); - ws.arr_fftw = fftw_alloc_real(static_cast<size_t>(ws.h * ws.w_fftw) * sizeof(fftw_complex)); + ws.arr_real = fftw_alloc_real(ws.h * ws.w_real); + ws.arr_fftw = fftw_alloc_real(ws.h * ws.w_fftw * sizeof(fftw_complex)); } void FourierTransform::init_r2c(int h, int w_real) @@ -281,7 +263,7 @@ void FourierTransform::rfft2complex_vec(complex2d_t& result) ASSERT(ws.arr_fftw); result.clear(); - result.resize(ws.h, complex1d_t(static_cast<size_t>(ws.w_fftw))); + result.resize(ws.h, complex1d_t(ws.w_fftw)); double* ptr = ws.arr_fftw; @@ -301,7 +283,7 @@ void FourierTransform::irfft2double_vec(double2d_t& result) result.clear(); // Resize the array for holding the FT output to correct dimensions - result.resize(ws.h, std::vector<double>(static_cast<size_t>(ws.w_real))); + result.resize(ws.h, std::vector<double>(ws.w_real)); double factor = ws.h * ws.w_real; @@ -318,7 +300,7 @@ void FourierTransform::fft2amp(complex2d_t& source, double2d_t& result) ASSERT(ws.arr_fftw); result.clear(); - result.resize(ws.h, std::vector<double>(static_cast<size_t>(ws.w_real))); + result.resize(ws.h, std::vector<double>(ws.w_real)); // Storing FT output for (int i = 0; i < ws.h; i++) { @@ -345,9 +327,9 @@ void FourierTransform::fftw_forward_FT(const double2d_t& src) *ptr = 0.0; // Building the input signal for fftw algorithm - for (size_t row = 0; row < static_cast<size_t>(ws.h); ++row) - for (size_t col = 0; col < static_cast<size_t>(ws.w_real); ++col) - ws.arr_real[static_cast<int>(row) * ws.w_real + static_cast<int>(col)] += src[row][col]; + for (int row = 0; row < ws.h; ++row) + for (int col = 0; col < ws.w_real; ++col) + ws.arr_real[row * ws.w_real + col] += src[row][col]; // Computing the FFT with fftw plan fftw_execute(ws.p_forw); @@ -362,8 +344,8 @@ void FourierTransform::fftw_backward_FT(const complex2d_t& src) // Building the fft input signal double* ptr = ws.arr_fftw; - for (size_t row = 0; row < static_cast<size_t>(ws.h); ++row) - for (size_t col = 0; col < static_cast<size_t>(ws.w_fftw); ++col) { + for (int row = 0; row < ws.h; ++row) + for (int col = 0; col < ws.w_fftw; ++col) { *ptr += src[row][col].real(); *(ptr + 1) += src[row][col].imag(); ptr += 2;