diff --git a/Device/Resolution/ConvolutionDetectorResolution.cpp b/Device/Resolution/ConvolutionDetectorResolution.cpp index 7ef0696ec31ad8964e71cd4d404fc7fca52acb0f..7b66916b7f6fd9864ffd0f9623550a68ff070cb7 100644 --- a/Device/Resolution/ConvolutionDetectorResolution.cpp +++ b/Device/Resolution/ConvolutionDetectorResolution.cpp @@ -101,55 +101,51 @@ void ConvolutionDetectorResolution::apply2dConvolution(Datafield* intensity_map) { ASSERT(m_res_function_2d); 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(); - size_t axis_size_2 = axis_2.size(); - if (axis_size_1 < 2 || axis_size_2 < 2) - return; // No 2d convolution for 1d data + const Scale& X = intensity_map->axis(0); + const Scale& Y = intensity_map->axis(1); + size_t nx = X.size(); + size_t ny = Y.size(); + ASSERT(nx > 1); + ASSERT(ny > 1); + // Construct source vector array from original 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(); - 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); - } + ASSERT(raw_data_size == nx * ny); + for (auto it = raw_source_vector.begin(); it != raw_source_vector.end(); it += ny) + source.emplace_back(std::vector<double>(it, it + ny)); + // Construct kernel vector from resolution function std::vector<std::vector<double>> kernel; - kernel.resize(axis_size_1); - double mid_value_1 = - axis_1.binCenter(axis_size_1 / 2); // because Convolve expects zero at midpoint - double mid_value_2 = - axis_2.binCenter(axis_size_2 / 2); // because Convolve expects zero at midpoint - double step_size_1 = - std::abs(axis_1.binCenter(0) - axis_1.binCenter(axis_size_1 - 1)) / (axis_size_1 - 1); - double step_size_2 = - std::abs(axis_2.binCenter(0) - axis_2.binCenter(axis_size_2 - 1)) / (axis_size_2 - 1); - for (size_t index_1 = 0; index_1 < axis_size_1; ++index_1) { - double value_1 = axis_1.binCenter(index_1) - mid_value_1; + kernel.resize(nx); + double mid_value1 = X.binCenter(nx / 2); // because Convolve expects zero at midpoint + double mid_value2 = Y.binCenter(ny / 2); // because Convolve expects zero at midpoint + double dx = std::abs(X.binCenter(0) - X.binCenter(nx - 1)) / (nx - 1); + double dy = std::abs(Y.binCenter(0) - Y.binCenter(ny - 1)) / (ny - 1); + for (size_t ix = 0; ix < nx; ++ix) { + double x = X.binCenter(ix) - mid_value1; std::vector<double> row_vector; - row_vector.resize(axis_size_2, 0.0); - for (size_t index_2 = 0; index_2 < axis_size_2; ++index_2) { - double value_2 = axis_2.binCenter(index_2) - mid_value_2; - double z_value = getIntegratedPDF2d(value_1, step_size_1, value_2, step_size_2); - row_vector[index_2] = z_value; + row_vector.resize(ny, 0.0); + for (size_t iy = 0; iy < ny; ++iy) { + double y = Y.binCenter(iy) - mid_value2; + row_vector[iy] = getIntegratedPDF2d(x, dx, y, dy); } - kernel[index_1] = row_vector; + kernel[ix] = row_vector; } // Calculate convolution std::vector<std::vector<double>> result; Convolve().fftconvolve(source, kernel, result); + // Populate intensity map with results std::vector<double> result_vector; - for (size_t index_1 = 0; index_1 < axis_size_1; ++index_1) { - for (size_t index_2 = 0; index_2 < axis_size_2; ++index_2) { - double value = result[index_1][index_2]; + for (size_t index1 = 0; index1 < nx; ++index1) { + for (size_t index2 = 0; index2 < ny; ++index2) { + double value = result[index1][index2]; result_vector.push_back(value); } } - ASSERT(axis_size_1 * axis_size_2 == intensity_map->size()); + ASSERT(nx * ny == intensity_map->size()); for (size_t i = 0; i < intensity_map->size(); ++i) { size_t i0 = intensity_map->frame().projectedIndex(i, 0); size_t i1 = intensity_map->frame().projectedIndex(i, 1);