diff --git a/Device/Resolution/ConvolutionDetectorResolution.cpp b/Device/Resolution/ConvolutionDetectorResolution.cpp index bfea137c8552b2bd22d6c2c25c6e33ef6c068ace..32614837e657772f4470f425629fbb6cac2a8094 100644 --- a/Device/Resolution/ConvolutionDetectorResolution.cpp +++ b/Device/Resolution/ConvolutionDetectorResolution.cpp @@ -17,7 +17,6 @@ #include "Base/Axis/Scale.h" #include "Base/Util/Assert.h" #include "Device/Resolution/Convolve.h" -#include <stdexcept> ConvolutionDetectorResolution::ConvolutionDetectorResolution(cumulative_DF_1d res_function_1d) : m_rank(1) @@ -56,23 +55,13 @@ std::vector<const INode*> ConvolutionDetectorResolution::nodeChildren() const void ConvolutionDetectorResolution::applyDetectorResolution(Datafield* intensity_map) const { - if (intensity_map->rank() != m_rank) { - throw std::runtime_error( - "ConvolutionDetectorResolution::applyDetectorResolution -> Error! " - "Intensity map must have same dimension as detector resolution function."); - } - switch (m_rank) { - case 1: + ASSERT(intensity_map->rank() == m_rank); + if (m_rank == 1) apply1dConvolution(intensity_map); - break; - case 2: + else if (m_rank == 2) apply2dConvolution(intensity_map); - break; - default: - throw std::runtime_error( - "ConvolutionDetectorResolution::applyDetectorResolution -> Error! " - "Class ConvolutionDetectorResolution must be initialized with dimension 1 or 2."); - } + else + ASSERT_NEVER; } void ConvolutionDetectorResolution::setResolutionFunction(const IResolutionFunction2D& resFunc) @@ -82,11 +71,9 @@ void ConvolutionDetectorResolution::setResolutionFunction(const IResolutionFunct void ConvolutionDetectorResolution::apply1dConvolution(Datafield* intensity_map) const { - ASSERT(m_res_function_1d == nullptr); - if (intensity_map->rank() != 1) - throw std::runtime_error( - "ConvolutionDetectorResolution::apply1dConvolution -> Error! " - "Number of axes for intensity map does not correspond to the dimension of the map."); + ASSERT(m_res_function_1d); + ASSERT(intensity_map->rank() == 1); + const Scale& axis = intensity_map->axis(0); // Construct source vector from original intensity map std::vector<double> source_vector = intensity_map->flatVector(); @@ -94,10 +81,7 @@ void ConvolutionDetectorResolution::apply1dConvolution(Datafield* intensity_map) if (data_size < 2) return; // No convolution for sets of zero or one element // Construct kernel vector from resolution function - if (axis.size() != data_size) - throw std::runtime_error( - "ConvolutionDetectorResolution::apply1dConvolution -> Error! " - "Size of axis for intensity map does not correspond to size of data in the map."); + ASSERT(axis.size() == data_size); double step_size = std::abs(axis.binCenter(0) - axis.binCenter(axis.size() - 1)) / (data_size - 1); double mid_value = axis.binCenter(axis.size() / 2); // because Convolve expects zero at midpoint @@ -117,10 +101,7 @@ void ConvolutionDetectorResolution::apply1dConvolution(Datafield* intensity_map) void ConvolutionDetectorResolution::apply2dConvolution(Datafield* intensity_map) const { ASSERT(m_res_function_2d); - if (intensity_map->rank() != 2) - throw std::runtime_error( - "ConvolutionDetectorResolution::apply2dConvolution -> Error! " - "Number of axes for intensity map does not correspond to the dimension of the map."); + ASSERT(intensity_map->rank() == 2); const Scale& axis_1 = intensity_map->axis(0); const Scale& axis_2 = intensity_map->axis(1); size_t axis_size_1 = axis_1.size(); @@ -131,10 +112,7 @@ void ConvolutionDetectorResolution::apply2dConvolution(Datafield* intensity_map) std::vector<double> raw_source_vector = intensity_map->flatVector(); std::vector<std::vector<double>> source; size_t raw_data_size = raw_source_vector.size(); - if (raw_data_size != axis_size_1 * axis_size_2) - throw std::runtime_error( - "ConvolutionDetectorResolution::apply2dConvolution -> Error! " - "Intensity map data size does not match the product of its axes' sizes"); + ASSERT(raw_data_size == axis_size_1 * axis_size_2); for (auto it = raw_source_vector.begin(); it != raw_source_vector.end(); it += axis_size_2) { std::vector<double> row_vector(it, it + axis_size_2); source.push_back(row_vector); diff --git a/Device/Resolution/Convolve.cpp b/Device/Resolution/Convolve.cpp index 5792a038e757321bcaaac5e393092a8214e8aed0..85c4de7dd779a64794a9dff372b6b0412e745642 100644 --- a/Device/Resolution/Convolve.cpp +++ b/Device/Resolution/Convolve.cpp @@ -13,9 +13,7 @@ // ************************************************************************************************ #include "Device/Resolution/Convolve.h" -#include <iostream> -#include <sstream> -#include <stdexcept> // need overlooked by g++ 5.4 +#include "Base/Util/Assert.h" Convolve::Convolve() : m_mode(FFTW_UNDEFINED) @@ -99,12 +97,11 @@ void Convolve::fftconvolve(const double2d_t& source, const double2d_t& kernel, d result[i].resize(ws.w_dst, 0); for (int j = 0; j < ws.w_dst; j++) { // result[i][j]=ws.dst[i*ws.w_dst+j]; - if (m_mode == FFTW_CIRCULAR_SAME_SHIFTED) { + if (m_mode == FFTW_CIRCULAR_SAME_SHIFTED) result[i][j] = ws.dst_fft[((i + int(ws.h_kernel / 2.0)) % ws.h_fftw) * ws.w_fftw + (j + int(ws.w_kernel / 2.0)) % ws.w_fftw]; - } else { + else result[i][j] = ws.dst_fft[(i + ws.h_offset) * ws.w_fftw + j + ws.w_offset]; - } } } } @@ -121,8 +118,7 @@ void Convolve::fftconvolve(const double1d_t& source, const double1d_t& kernel, d double2d_t result2d; fftconvolve(source2d, kernel2d, result2d); - if (result2d.size() != 1) - throw std::runtime_error("Convolve::fftconvolve -> Panic in 1d"); + ASSERT(result2d.size() == 1); result = result2d[0]; } @@ -131,12 +127,10 @@ void Convolve::fftconvolve(const double1d_t& source, const double1d_t& kernel, d /* ************************************************************************* */ void Convolve::init(int h_src, int w_src, int h_kernel, int w_kernel) { - if (!h_src || !w_src || !h_kernel || !w_kernel) { - std::ostringstream os; - os << "Convolve::init -> Panic! Wrong dimensions " << h_src << " " << w_src << " " - << h_kernel << " " << w_kernel << std::endl; - throw std::runtime_error(os.str()); - } + ASSERT(h_src); + ASSERT(w_src); + ASSERT(h_kernel); + ASSERT(w_kernel); ws.clear(); ws.h_src = h_src; @@ -178,19 +172,12 @@ void Convolve::init(int h_src, int w_src, int h_kernel, int w_kernel) break; case FFTW_LINEAR_VALID: // Valid Linear convolution - if (ws.h_kernel > ws.h_src || ws.w_kernel > ws.w_src) { - ws.h_fftw = 0; - ws.w_fftw = 0; - ws.h_dst = 0; - ws.w_dst = 0; - std::cout << "The 'valid' convolution results in an empty matrix" << std::endl; - throw std::runtime_error("The 'valid' convolution results in an empty matrix"); - } else { - ws.h_fftw = find_closest_factor(h_src); - ws.w_fftw = find_closest_factor(w_src); - ws.h_dst = h_src - h_kernel + 1; - ws.w_dst = w_src - w_kernel + 1; - } + ASSERT(ws.h_kernel <= ws.h_src); + ASSERT(ws.w_kernel <= ws.w_src); + ws.h_fftw = find_closest_factor(h_src); + ws.w_fftw = find_closest_factor(w_src); + ws.h_dst = h_src - h_kernel + 1; + ws.w_dst = w_src - w_kernel + 1; // Here we just take [h_dst x w_dst] elements starting at [h_kernel-1;w_kernel-1] ws.h_offset = ws.h_kernel - 1; ws.w_offset = ws.w_kernel - 1; @@ -210,11 +197,7 @@ void Convolve::init(int h_src, int w_src, int h_kernel, int w_kernel) ws.w_offset = 0; break; default: - std::cout - << "Unrecognized convolution mode, possible modes are " - << "FFTW_LINEAR_FULL, FFTW_LINEAR_SAME, FFTW_LINEAR_SAME_UNPADDED, FFTW_LINEAR_VALID, " - << "FFTW_CIRCULAR_SAME, FFTW_CIRCULAR_SHIFTED " << std::endl; - break; + ASSERT_NEVER; } ws.in_src = new double[ws.h_fftw * ws.w_fftw]; @@ -228,19 +211,15 @@ void Convolve::init(int h_src, int w_src, int h_kernel, int w_kernel) // Initialization of the plans ws.p_forw_src = fftw_plan_dft_r2c_2d(ws.h_fftw, ws.w_fftw, ws.in_src, (fftw_complex*)ws.out_src, FFTW_ESTIMATE); - if (ws.p_forw_src == nullptr) - throw std::runtime_error("Convolve::init -> Error! Cannot initialise p_forw_src plan."); - + ASSERT(ws.p_forw_src); ws.p_forw_kernel = fftw_plan_dft_r2c_2d(ws.h_fftw, ws.w_fftw, ws.in_kernel, (fftw_complex*)ws.out_kernel, FFTW_ESTIMATE); - if (ws.p_forw_kernel == nullptr) - throw std::runtime_error("Convolve::init -> Error! Cannot initialise p_forw_kernel plan."); + ASSERT(ws.p_forw_kernel); // The backward FFT takes ws.out_kernel as input ws.p_back = fftw_plan_dft_c2r_2d(ws.h_fftw, ws.w_fftw, (fftw_complex*)ws.out_kernel, ws.dst_fft, FFTW_ESTIMATE); - if (ws.p_back == nullptr) - throw std::runtime_error("Convolve::init -> Error! Cannot initialise p_back plan."); + ASSERT(ws.p_back); } /* ************************************************************************* */ @@ -249,8 +228,8 @@ void Convolve::init(int h_src, int w_src, int h_kernel, int w_kernel) void Convolve::fftw_circular_convolution(const double2d_t& src, const double2d_t& kernel) { - if (ws.h_fftw <= 0 || ws.w_fftw <= 0) - throw std::runtime_error("Convolve::fftw_convolve -> Panic! Initialization is missed."); + ASSERT(ws.h_fftw > 0); + ASSERT(ws.w_fftw > 0); double *ptr, *ptr_end, *ptr2; @@ -319,9 +298,8 @@ bool Convolve::is_optimal(int n) if (n == 1) return false; size_t ntest = n; - for (size_t i = 0; i < m_implemented_factors.size(); i++) { + for (size_t i = 0; i < m_implemented_factors.size(); i++) while ((ntest % m_implemented_factors[i]) == 0) ntest = ntest / m_implemented_factors[i]; - } return ntest == 1; }