3
\$\begingroup\$

This is a follow-up question for im2double and im2uint8 Functions Implementation for Image in C++, conv2 Template Function Implementation for Image in C++ and An Updated Multi-dimensional Image Data Structure with Variadic Template Functions in C++. I am trying to implement imgaussfilt template function like Matlab's imgaussfilt to perform Gaussian blur in this post.

The experimental implementation

  • imgaussfilt template function implementation (in file image_operations.h)

    namespace TinyDIP
{
    //  imgaussfilt template function implementation
    template<typename ElementT>
    constexpr static auto imgaussfilt(const Image<ElementT>& input, double sigma = 0.5, bool is_size_same = true)
    {
        if (input.getDimensionality()!=2)
        {
            throw std::runtime_error("Unsupported dimension!");
        }
        return imgaussfilt(input, sigma, static_cast<int>(2 * std::ceil(2 * sigma) + 1), is_size_same);
    }

    //  imgaussfilt template function implementation
    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, bool is_size_same = true)
    {
        if (input.getDimensionality()!=2)
        {
            throw std::runtime_error("Unsupported dimension!");
        }
        return imgaussfilt(input, sigma, sigma, filter_size, is_size_same);
    }

    //  imgaussfilt template function implementation
    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 sigma1, SigmaT sigma2, SizeT filter_size, bool is_size_same = true)
    {
        if (input.getDimensionality()!=2)
        {
            throw std::runtime_error("Unsupported dimension!");
        }
        auto filter_mask = gaussianFigure2D(
                                    filter_size,
                                    filter_size,
                                    (static_cast<double>(filter_size) + 1.0) / 2.0,
                                    (static_cast<double>(filter_size) + 1.0) / 2.0,
                                    sigma1, sigma2);
        auto sum_result = sum(filter_mask);
        filter_mask = divides(filter_mask, sum_result);             //  Normalization
        return conv2(input, filter_mask, is_size_same);
    }
}

  • gaussianFigure2D template function implementation (in file image_operations.h)

    namespace TinyDIP
{
    //  multiple standard deviations
    template<class InputT>
    constexpr static Image<InputT> gaussianFigure2D(
        const std::size_t xsize, const std::size_t ysize, 
        const std::size_t centerx, const std::size_t centery,
        const InputT standard_deviation_x, const InputT standard_deviation_y)
    {
        auto output = Image<InputT>(xsize, ysize);
        auto row_vector_x = Image<InputT>(xsize, std::size_t{1});
        for (std::size_t x = 0; x < xsize; ++x)
        {
            row_vector_x.at(x, 0) = normalDistribution1D(static_cast<InputT>(x) - static_cast<InputT>(centerx), standard_deviation_x);
        }

        auto row_vector_y = Image<InputT>(ysize, std::size_t{1});
        for (std::size_t y = 0; y < ysize; ++y)
        {
            row_vector_y.at(y, 0) = normalDistribution1D(static_cast<InputT>(y) - static_cast<InputT>(centery), standard_deviation_y);
        }

        for (std::size_t y = 0; y < ysize; ++y)
        {
            for (std::size_t x = 0; x < xsize; ++x)
            {
                output.at(x, y) = row_vector_x.at(x, 0) * row_vector_y.at(y, 0);
            }
        }
        return output;
    }
}

  • conv2 template function implementation (in file image_operations.h)

    namespace TinyDIP
{
    //  conv2 template function implementation
    template<typename ElementT>
    requires(std::floating_point<ElementT> || std::integral<ElementT> || is_complex<ElementT>::value)
    constexpr auto conv2(const Image<ElementT>& x, const Image<ElementT>& y, bool is_size_same = false)
    {
        auto output = Image<ElementT>(x.getWidth() + y.getWidth() - 1, x.getHeight() + y.getHeight() - 1);
        for (std::size_t y1 = 0; y1 < x.getHeight(); ++y1) {
            auto* x_row = &(x.at(0, y1));
            for (std::size_t y2 = 0; y2 < y.getHeight(); ++y2) {
                auto* y_row = &(y.at(0, y2));
                auto* out_row = &(output.at(0, y1 + y2));
                for (std::size_t x1 = 0; x1 < x.getWidth(); ++x1) {
                    for (std::size_t x2 = 0; x2 < y.getWidth(); ++x2) {
                        out_row[x1 + x2] += x_row[x1] * y_row[x2];
                    }
                }
            }
        }
        if(is_size_same)
        {
            output = subimage(output, x.getWidth(), x.getHeight(), static_cast<double>(output.getWidth()) / 2.0, static_cast<double>(output.getHeight()) / 2.0);
        }
        return output;
    }

    //  conv2 template function implementation
    template<typename ElementT, typename ElementT2>
    requires (((std::same_as<ElementT, RGB>) || (std::same_as<ElementT, RGB_DOUBLE>) || (std::same_as<ElementT, HSV>)) &&
              (std::floating_point<ElementT2> || std::integral<ElementT2> || is_complex<ElementT2>::value))
    constexpr static auto conv2(const Image<ElementT>& input1, const Image<ElementT2>& input2, bool is_size_same = false)
    {
        return apply_each(input1, [&](auto&& planes) { return conv2(planes, input2, is_size_same); });
    }
}

The usage of imgaussfilt function:

void imgaussfiltTest(std::string_view input_image_path = "InputImages/1", std::string_view output_image_path = "OutputImages/imgaussfiltTest")
{
    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::imgaussfilt(TinyDIP::im2double(input_img), sigma)
                                );
        TinyDIP::bmp_write(
            (std::string(output_image_path) + std::string("_sigma=") + std::to_string(sigma)).c_str(),
            output_img);
    }
    
}

int main(int argc, char* argv[])
{
    auto start = std::chrono::system_clock::now();
    imgaussfiltTest();
    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;
}

TinyDIP on GitHub

All suggestions are welcome.

The summary information:

\$\endgroup\$

1 Answer 1

Reset to default
4
\$\begingroup\$

I have strong opinions about the Gaussian filter, don't take this personally. :)

First of all, you should know that the Gaussian filter is separable. That is, you can apply a 1D Gaussian filter to the rows of an image, and then another 1D Gaussian filter to the columns, to obtain a 2D Gaussian convolution. There are several advantages to implementing the filter this way:

  • It is significantly cheaper computationally (2*filter_size operations per pixel instead of filter_size*filter_size operations per pixel).
  • It is trivial to implement the filter for an image of arbitrary dimensionality.
  • It becomes fairly cheap to do more meaningful padding than the implicit zero-padding you're using in conv2().

I've written some blog post in the past about Gaussian filtering, please read them:

API

You have three signatures:

  • imgaussfilt(const Image<ElementT>& input, double sigma = 0.5, bool is_size_same = true)
  • imgaussfilt(const Image<ElementT>& input, SigmaT sigma, SizeT filter_size, bool is_size_same = true)
  • imgaussfilt(const Image<ElementT>& input, SigmaT sigma1, SigmaT sigma2, SizeT filter_size, bool is_size_same = true)

These exist only to avoid input arguments with default values. For example, the first two could be merged by giving filter_size a default value of 0, and having the function compute an appropriate size unless the user specifies a positive value.

The third one is dangerous: the 2D kernel is always square, but the two sigmas can be different. You're relying on the user knowing how to set this size properly. And a correct size value here would depend on the larger of the two sigmas, leading to a kernel potentially with many zeros in it. For example, say I set the sigmas to 5 and 1, the kernel will be 5 times larger than necessary in one dimension.

There is also no signature where I can give two different sigmas, and have the library compute the right kernel sizes for me.

But in general, I would prefer not to offer the filter_size parameter at all, as it is easy to misuse due to lack of knowledge. See for example the third blog post I linked above, where a scientific paper from a major university used "an 11x11 pixel Gaussian filter with standard deviation σ=25 pixels", which is not a Gaussian. Having a truncation parameter instead makes the function much safer. Having the user set truncation=0.3 requires them to really know that they're doing something weird. Your function would then compute filter_size = static_cast<int>(2 * std::ceil(truncation * sigma) + 1).

Also, in my world, truncation defaults to 3, 2 is typically not enough. See again that third blog post above. I just learned that imgaussfilt in MATLAB uses 2, I think that's an unfortunate choice.

As someone that uses Gaussian filters a lot, I would prefer a function with a signature like this:

imgaussfilt(const Image<ElementT>& input, std::vector<SigmaT> sigmas, BoundaryCondition boundaryCondition = BoundaryCondition::mirror, double truncation = 3)

sigmas is an array, so that I can specify one value per image dimension. If given a single value, that value would be used for all dimensions.

The BoundaryCondition enum would contain the various possible ways to handle pixels outside the image domain. This is an important feature to be able to control in specialized circumstances. In contrast, the is_size_same option I would have no use for.

Code issues

What happens if somebody sets filter_size = 0?

Possible extensions

This filter always uses zero padding as a boundary condition (conv2 implicitly assumes zero values outside the image domain). This is often a poor choice (though not always, sometimes it is exactly what you need). The mirror boundary condition (sometimes called reflect) tends to produce more useful results near the image boundary. To implement this in conv2 you'd have to copy the image into a larger array, and fill the new areas by copy. But in a separable convolution, where the operation is applied to 1D lines across the image, it is quite cheap to copy one image line into a buffer, extend this line on both sides using the chosen boundary condition, and apply the 1D filter to that. One image line plus its padding usually will fit in the cache, making the intermediate line buffer not too expensive.

Another advantage of the intermediate line buffer is that you can then write the result of the filter back into the input image. This makes the repeated 1D filtering used by the separable filter much more memory friendly.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.