This is a follow-up question for conv2 Template Function Implementation for Image in C++ and imgaussfilt Template Function Implementation for Image in C++. I am trying to perform Difference of Gaussians with imgaussfilt template function in this post. Here, imgaussfilt template function is enhanced with BoundaryCondition option, which can be set to constant (pad image with elements of constant value), mirror (pad image with mirror reflections of itself) or replicate (pad by repeating border elements of array). An example image output from difference_of_gaussian template function:
The experimental implementation
difference_of_gaussiantemplate function implementation (in fileimage_operations.h)namespace TinyDIP { // difference_of_gaussian template function implementation template<typename ElementT, typename SigmaT = double> requires(std::floating_point<SigmaT> || std::integral<SigmaT>) constexpr static auto difference_of_gaussian( const Image<ElementT>& input, SigmaT sigma1, SigmaT sigma2) { return subtract( imgaussfilt(input, sigma1, static_cast<int>(computeFilterSizeFromSigma(sigma1))), imgaussfilt(input, sigma2, static_cast<int>(computeFilterSizeFromSigma(sigma2))) ); } }computeFilterSizeFromSigmatemplate function implementation (in fileimage_operations.h)namespace TinyDIP { // computeFilterSizeFromSigma template function implementation template<typename ElementT> constexpr static auto computeFilterSizeFromSigma(ElementT sigma) { return 2 * std::ceil(2 * sigma) + 1; } }imgaussfilttemplate function implementation (in fileimage_operations.h)namespace TinyDIP { // imgaussfilt template function implementation // https://codereview.stackexchange.com/q/292985/231235 // giving filter_size a default value of 0, and having the function compute an appropriate size unless the user specifies a positive value. template<typename ElementT, typename SigmaT = double, std::integral SizeT = int> requires(std::floating_point<SigmaT> || std::integral<SigmaT>) constexpr static auto imgaussfilt( const Image<ElementT>& input, SigmaT sigma, SizeT filter_size = 0, BoundaryCondition boundaryCondition = BoundaryCondition::mirror, ElementT value_for_constant_padding = ElementT{}) { if (input.getDimensionality()!=2) { throw std::runtime_error("Unsupported dimension!"); } if (filter_size == 0) { return imgaussfilt( std::execution::seq, input, sigma, sigma, static_cast<int>(computeFilterSizeFromSigma(sigma)), static_cast<int>(computeFilterSizeFromSigma(sigma)), boundaryCondition, value_for_constant_padding); } return imgaussfilt( std::execution::seq, input, sigma, sigma, filter_size, filter_size, boundaryCondition, value_for_constant_padding); } // imgaussfilt template function implementation // https://codereview.stackexchange.com/q/292985/231235 template<class ExecutionPolicy, typename ElementT, typename SigmaT = double, std::integral SizeT = int> requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>)&& (std::floating_point<SigmaT> || std::integral<SigmaT>) constexpr static auto imgaussfilt( ExecutionPolicy&& execution_policy, const Image<ElementT>& input, SigmaT sigma1, SigmaT sigma2, SizeT filter_size1, SizeT filter_size2, BoundaryCondition boundaryCondition = BoundaryCondition::mirror, ElementT value_for_constant_padding = ElementT{}) { if (input.getDimensionality()!=2) { throw std::runtime_error("Unsupported dimension!"); } Image<ElementT> padded_image; switch(boundaryCondition) { case constant: padded_image = generate_constant_padding_image(execution_policy, input, filter_size1, filter_size2, value_for_constant_padding); break; case mirror: padded_image = generate_mirror_padding_image(execution_policy, input, filter_size1, filter_size2, value_for_constant_padding); break; case replicate: padded_image = generate_replicate_padding_image(execution_policy, input, filter_size1, filter_size2, value_for_constant_padding); break; } auto filter_mask_x = gaussianFigure1D( filter_size1, (static_cast<double>(filter_size1) + 1.0) / 2.0, sigma1); auto sum_result = sum(filter_mask_x); filter_mask_x = divides(filter_mask_x, sum_result); // Normalization auto output = conv2(padded_image, filter_mask_x, true); auto filter_mask_y = transpose(gaussianFigure1D( filter_size2, (static_cast<double>(filter_size2) + 1.0) / 2.0, sigma2)); sum_result = sum(filter_mask_y); filter_mask_y = divides(filter_mask_y, sum_result); // Normalization output = conv2(output, filter_mask_y, true); output = subimage(output, input.getWidth(), input.getHeight(), static_cast<double>(output.getWidth()) / 2.0, static_cast<double>(output.getHeight()) / 2.0); return output; } }BoundaryConditionenumeration declarationenum BoundaryCondition { constant, mirror, replicate };generate_constant_padding_imagetemplate function implementation (in fileimage_operations.h)namespace TinyDIP { // generate_constant_padding_image template function implementation template<typename ElementT> constexpr static auto generate_constant_padding_image( const Image<ElementT> input, std::size_t width_expansion, std::size_t height_expansion, ElementT default_value = ElementT{}) { return generate_constant_padding_image(std::execution::seq, input, width_expansion, height_expansion, default_value); } // generate_constant_padding_image template function implementation (with Execution Policy) template<class ExecutionPolicy, typename ElementT> requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>) constexpr static auto generate_constant_padding_image( ExecutionPolicy&& execution_policy, const Image<ElementT> input, std::size_t width_expansion, std::size_t height_expansion, ElementT default_value = ElementT{}) { if (input.getDimensionality()!=2) { throw std::runtime_error("Unsupported dimension!"); } Image<ElementT> output(input.getWidth() + 2 * width_expansion, input.getHeight() + 2 * height_expansion); output.setAllValue(default_value); output = paste2D(execution_policy, output, input, width_expansion, height_expansion, default_value); return output; } }generate_mirror_padding_imagetemplate function implementation (in fileimage_operations.h)namespace TinyDIP { // generate_mirror_padding_image template function implementation template<typename ElementT> constexpr static auto generate_mirror_padding_image( const Image<ElementT> input, std::size_t width_expansion, std::size_t height_expansion, ElementT default_value = ElementT{}) { if (input.getDimensionality()!=2) { throw std::runtime_error("Unsupported dimension!"); } return generate_mirror_padding_image(std::execution::seq, input, width_expansion, height_expansion, default_value); } // generate_mirror_padding_image template function implementation (with Execution Policy) template<class ExecutionPolicy, typename ElementT> requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>) constexpr static auto generate_mirror_padding_image( ExecutionPolicy&& execution_policy, const Image<ElementT> input, std::size_t width_expansion, std::size_t height_expansion, ElementT default_value = ElementT{}) { if (input.getDimensionality()!=2) { throw std::runtime_error("Unsupported dimension!"); } auto output = generate_constant_padding_image(execution_policy, input, width_expansion, height_expansion, default_value); auto flipped_vertical = flip_vertical(input); output = paste2D( execution_policy, output, subimage2(flipped_vertical, 0, flipped_vertical.getWidth() - 1, input.getHeight() - height_expansion - 1, flipped_vertical.getHeight() - 1), width_expansion, 0, default_value); output = paste2D( execution_policy, output, subimage2(flipped_vertical, 0, flipped_vertical.getWidth() - 1, 0, height_expansion), width_expansion, input.getHeight() + height_expansion - 1, default_value); auto flipped_horizontal = flip_horizontal(input); output = paste2D( execution_policy, output, subimage2(flipped_horizontal, input.getWidth() - width_expansion - 1, flipped_horizontal.getWidth() - 1, 0, flipped_horizontal.getHeight() - 1), 0, height_expansion, default_value); output = paste2D( execution_policy, output, subimage2(flipped_horizontal, 0, width_expansion, 0, flipped_horizontal.getHeight() - 1), input.getWidth() + width_expansion - 1, height_expansion, default_value); auto flipped_horizontal_vertical = flip_horizontal_vertical(input); output = paste2D( execution_policy, output, subimage2( flipped_horizontal_vertical, flipped_horizontal_vertical.getWidth() - width_expansion - 1, flipped_horizontal_vertical.getWidth() - 1, flipped_horizontal_vertical.getHeight() - height_expansion - 1, flipped_horizontal_vertical.getHeight() - 1), 0, 0, default_value); output = paste2D( execution_policy, output, subimage2( flipped_horizontal_vertical, 0, width_expansion, flipped_horizontal_vertical.getHeight() - height_expansion - 1, flipped_horizontal_vertical.getHeight() - 1), input.getWidth() + width_expansion - 1, 0, default_value); output = paste2D( execution_policy, output, subimage2( flipped_horizontal_vertical, flipped_horizontal_vertical.getWidth() - width_expansion - 1, flipped_horizontal_vertical.getWidth() - 1, 0, height_expansion), 0, input.getHeight() + height_expansion - 1, default_value); output = paste2D( execution_policy, output, subimage2( flipped_horizontal_vertical, 0, width_expansion, 0, height_expansion), input.getWidth() + width_expansion - 1, input.getHeight() + height_expansion - 1, default_value); return output; } }generate_replicate_padding_imagetemplate function implementation (in fileimage_operations.h)namespace TinyDIP { // generate_replicate_padding_image template function implementation template<typename ElementT> constexpr static auto generate_replicate_padding_image( const Image<ElementT> input, std::size_t width_expansion, std::size_t height_expansion, ElementT default_value = ElementT{}) { if (input.getDimensionality()!=2) { throw std::runtime_error("Unsupported dimension!"); } return generate_replicate_padding_image(std::execution::seq, input, width_expansion, height_expansion, default_value); } // generate_replicate_padding_image template function implementation (with Execution Policy) // Test: https://godbolt.org/z/1hebz7hEh template<class ExecutionPolicy, typename ElementT> requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>) constexpr static auto generate_replicate_padding_image( ExecutionPolicy&& execution_policy, const Image<ElementT> input, std::size_t width_expansion, std::size_t height_expansion, ElementT default_value = ElementT{}) { if (input.getDimensionality()!=2) { throw std::runtime_error("Unsupported dimension!"); } auto output = generate_constant_padding_image(execution_policy, input, width_expansion, height_expansion, default_value); // Top block for(std::size_t y = 0; y < height_expansion; ++y) { output = paste2D( execution_policy, output, subimage2(input, 0, input.getWidth() - 1, 0, 0), width_expansion, y, default_value); } // Bottom block for(std::size_t y = input.getHeight() + height_expansion; y < input.getHeight() + 2 * height_expansion; ++y) { output = paste2D( execution_policy, output, subimage2(input, 0, input.getWidth() - 1, input.getHeight() - 1, input.getHeight() - 1), width_expansion, y, default_value); } // Left block for(std::size_t x = 0; x < width_expansion; ++x) { output = paste2D( execution_policy, output, subimage2(input, 0, 0, 0, input.getHeight() - 1), x, height_expansion, default_value); } // Right block for(std::size_t x = input.getWidth() + width_expansion; x < input.getWidth() + 2 * width_expansion; ++x) { output = paste2D( execution_policy, output, subimage2(input, input.getWidth() - 1, input.getWidth() - 1, 0, input.getHeight() - 1), x, height_expansion, default_value); } Image<ElementT> temp(width_expansion, height_expansion); // Left-top corner temp.setAllValue(input.at(0, 0)); output = paste2D( execution_policy, output, temp, 0, 0, default_value); // Right-top corner temp.setAllValue(input.at(input.getWidth() - 1, 0)); output = paste2D( execution_policy, output, temp, width_expansion + input.getWidth(), 0, default_value); // Left-bottom corner temp.setAllValue(input.at(0, input.getHeight() - 1)); output = paste2D( execution_policy, output, temp, 0, height_expansion + input.getHeight(), default_value); // Right-bottom corner temp.setAllValue(input.at(input.getWidth() - 1, input.getHeight() - 1)); output = paste2D( execution_policy, output, temp, width_expansion + input.getWidth(), height_expansion + input.getHeight(), default_value); return output; } }subtracttemplate function implementation (in fileimage_operations.h)namespace TinyDIP { // subtract Template Function Implementation template<class InputT> constexpr static Image<InputT> subtract(const Image<InputT>& input1, const Image<InputT>& input2) { check_size_same(input1, input2); return pixelwiseOperation(std::minus<>{}, input1, input2); } // subtract Template Function Implementation template<class InputT> constexpr static auto subtract(const std::vector<Image<InputT>>& input1, const std::vector<Image<InputT>>& input2) { assert(input1.size() == input2.size()); return recursive_transform<1>( [](auto&& input1_element, auto&& input2_element) { return subtract(input1_element, input2_element); }, input1, input2); } // subtract Function Implementation constexpr static Image<RGB> subtract(const Image<RGB>& input1, const Image<RGB>& input2) { check_size_same(input1, input2); return pixelwiseOperation( [](RGB x, RGB y) { RGB rgb; for(std::size_t channel_index = 0; channel_index < 3; ++channel_index) { rgb.channels[channel_index] = std::clamp( x.channels[channel_index] - y.channels[channel_index], 0, 255); } return rgb; }, input1, input2 ); } // subtract Template Function Implementation template<class InputT> requires((std::same_as<InputT, RGB_DOUBLE>) || (std::same_as<InputT, HSV>)) constexpr static auto subtract(const Image<InputT>& input1, const Image<InputT>& input2) { check_size_same(input1, input2); return pixelwiseOperation( [](InputT x, InputT y) { InputT output; for(std::size_t channel_index = 0; channel_index < 3; ++channel_index) { output.channels[channel_index] = x.channels[channel_index] - y.channels[channel_index]; } return output; }, input1, input2 ); } }paste2Dtemplate function implementation (in fileimage_operations.h)namespace TinyDIP { // paste2D template function implementation template<typename ElementT> constexpr static auto paste2D(const Image<ElementT>& background, const Image<ElementT>& target, std::size_t x_location, std::size_t y_location, ElementT default_value = ElementT{}) { return paste2D(std::execution::seq, background, target, x_location, y_location, default_value); } // paste2D template function implementation (with execution policy) // Test: https://godbolt.org/z/5hjns1nGP template<class ExecutionPolicy, typename ElementT> requires(std::is_execution_policy_v<std::remove_cvref_t<ExecutionPolicy>>) constexpr static auto paste2D(ExecutionPolicy&& execution_policy, const Image<ElementT>& background, const Image<ElementT>& target, std::size_t x_location, std::size_t y_location, ElementT default_value = ElementT{}) { if (background.getDimensionality()!=2) { throw std::runtime_error("Unsupported dimension!"); } if (target.getDimensionality()!=2) { throw std::runtime_error("Unsupported dimension!"); } if((background.getWidth() >= target.getWidth() + x_location) && (background.getHeight() >= target.getHeight() + y_location)) { auto output = background; for (std::size_t y = 0; y < target.getHeight(); ++y) { for (std::size_t x = 0; x < target.getWidth(); ++x) { output.at_without_boundary_check(x_location + x, y_location + y) = target.at_without_boundary_check(x, y); } } return output; } else { std::vector<ElementT> data; auto xsize = (background.getWidth() >= target.getWidth() + x_location)? background.getWidth(): (target.getWidth() + x_location); auto ysize = (background.getHeight() >= target.getHeight() + y_location)? background.getHeight(): (target.getHeight() + y_location); data.resize(xsize * ysize); std::fill(execution_policy, std::ranges::begin(data), std::ranges::end(data), default_value); Image<ElementT> output(data, xsize, ysize); for (std::size_t y = 0; y < background.getHeight(); ++y) { for (std::size_t x = 0; x < background.getWidth(); ++x) { output.at_without_boundary_check(x, y) = background.at_without_boundary_check(x, y); } } for (std::size_t y = 0; y < target.getHeight(); ++y) { for (std::size_t x = 0; x < target.getWidth(); ++x) { output.at_without_boundary_check(x_location + x, y_location + y) = target.at_without_boundary_check(x, y); } } return output; } } }subimage2template function implementation (in fileimage_operations.h)namespace TinyDIP { // subimage2 template function implementation template<typename ElementT> constexpr static auto subimage2(const Image<ElementT>& input, const std::size_t startx, const std::size_t endx, const std::size_t starty, const std::size_t endy) { assert(startx <= endx); assert(starty <= endy); Image<ElementT> output(endx - startx + 1, endy - starty + 1); for (std::size_t y = 0; y < output.getHeight(); ++y) { for (std::size_t x = 0; x < output.getWidth(); ++x) { output.at_without_boundary_check(x, y) = input.at_without_boundary_check(startx + x, starty + y); } } return output; } }flip_verticaltemplate function implementation (in fileimage_operations.h)namespace TinyDIP { // flip_vertical template function implementation template<typename ElementT> constexpr static auto flip_vertical(const Image<ElementT>& input) { if (input.getDimensionality()!=2) { throw std::runtime_error("Unsupported dimension!"); } Image<ElementT> output = input; #pragma omp parallel for collapse(2) for(std::size_t y = 0; y < input.getHeight(); ++y) { for(std::size_t x = 0; x < input.getWidth(); ++x) { output.at_without_boundary_check(x, input.getHeight() - y - 1) = input.at_without_boundary_check(x, y); } } return output; } }flip_horizontaltemplate function implementation (in fileimage_operations.h)namespace TinyDIP { // flip_horizontal template function implementation template<typename ElementT> constexpr static auto flip_horizontal(const Image<ElementT>& input) { if (input.getDimensionality()!=2) { throw std::runtime_error("Unsupported dimension!"); } Image<ElementT> output = input; #pragma omp parallel for collapse(2) for(std::size_t y = 0; y < input.getHeight(); ++y) { for(std::size_t x = 0; x < input.getWidth(); ++x) { output.at_without_boundary_check(input.getWidth() - x - 1, y) = input.at_without_boundary_check(x, y); } } return output; } }flip_horizontal_verticaltemplate function implementation (in fileimage_operations.h)namespace TinyDIP { // flip_horizontal_vertical template function implementation template<typename ElementT> constexpr static auto flip_horizontal_vertical(const Image<ElementT>& input) { if (input.getDimensionality()!=2) { throw std::runtime_error("Unsupported dimension!"); } Image<ElementT> output = input; #pragma omp parallel for collapse(2) for(std::size_t y = 0; y < input.getHeight(); ++y) { for(std::size_t x = 0; x < input.getWidth(); ++x) { output.at_without_boundary_check( input.getWidth() - x - 1, input.getHeight() - y - 1 ) = input.at_without_boundary_check(x, y); } } return output; } }
The usage of difference_of_gaussian template function:
void differenceOfGaussianTest(std::string_view input_image_path = "InputImages/1", std::string_view output_image_path = "OutputImages/differenceOfGaussianTest")
{
auto input_img = TinyDIP::bmp_read(std::string(input_image_path).c_str(), false);
for(int sigma = 1; sigma < 10; ++sigma)
{
auto output_img = TinyDIP::im2uint8(
TinyDIP::multiplies(
TinyDIP::abs(
TinyDIP::difference_of_gaussian(TinyDIP::im2double(input_img), static_cast<double>(sigma), static_cast<double>(sigma) - 1.0)
),
3
)
);
TinyDIP::bmp_write(
(std::string(output_image_path) + std::string("_sigma=") + std::to_string(sigma)).c_str(),
output_img);
}
return;
}
int main(int argc, char* argv[])
{
auto start = std::chrono::system_clock::now();
differenceOfGaussianTest();
auto end = std::chrono::system_clock::now();
std::chrono::duration<double> elapsed_seconds = end - start;
std::time_t end_time = std::chrono::system_clock::to_time_t(end);
std::cout << "Computation finished at " << std::ctime(&end_time) << "elapsed time: " << elapsed_seconds.count() << '\n';
return EXIT_SUCCESS;
}
All suggestions are welcome.
The summary information:
Which question it is a follow-up to?
conv2 Template Function Implementation for Image in C++ and
imgaussfilt Template Function Implementation for Image in C++
What changes has been made in the code since last question?
I am trying to perform Difference of Gaussians with
imgaussfilttemplate function in this post. Inimgaussfilttemplate function, padding options are added, and two 1D Gaussian filters are applied to rows and columns respectively of an image instead of 2D Gaussian filter to reduce the computation complexity.Why a new review is being asked for?
Please review the implementation of
difference_of_gaussiantemplate function and its tests.
subtract()is not defined for gray-scale images? Note that the difference of Gaussians has a signed result, your example output doesn't have any negative values, which means half of the pixels have their information missing. \$\endgroup\$paste2D()andsubimage2()functions, as well as theflip_vertical()andflip_horizontal()functions. They're at the core of the padding functions, which is the most messy thing to implement in here. \$\endgroup\$