1
votes

I'm trying to blur a matrix using Gaussian 2D Convolution. But I'm getting sharp transitions at boundary elements.

Here is a piece of code that I'm running:

// create 1D Kernel
void createGaussianKerenel_1D() {
    unsigned kernelSize = 2 * kernelRad_ + 1;
    gaussian1Dkernel_ = vector<double>(kernelSize);
    double sigma = (double)kernelRad_;

    double sum = 0.0;
    for(unsigned i = 0; i < kernelSize; ++i) {
        gaussian1Dkernel_[i] = gaussian(i, sigma);
        sum += gaussian1Dkernel_[i];
    }

    // normalize
    for(unsigned i = 0; i < kernelSize; ++i) {
        gaussian1Dkernel_[i] /= sum;
        cout << gaussian1Dkernel_[i] << endl;
    }
}

// gaussian function
double gaussian(unsigned int i, double sigma) const {
    double x = ((double)i - (double)kernelRad_) / sigma;

    return exp(-x * x / 2);
}

// do Separable 2D Convolution (in place)
// my initialMatrix_ is of yn_ x xn_ size
void getBlurredThermalMap() {
    assert(!gaussian1Dkernel_.empty());
    vector<vector<double> > tmpMatrix(yn_);
    unsigned kernelSize = 2 * kernelRad_ + 1;

    // in x direction
    for(unsigned i = 0; i < yn_; ++i) {
        for(unsigned j = 0; j < xn_; ++j) {
            double approxVal = 0.0;
            for(unsigned row = 0; row < kernelSize; ++row) {
                unsigned neighbor_j = j + row - kernelRad_;
                // ignore values that are out of bound
                if(neighbor_j >= 0 && neighbor_j < xn_) {
                    approxVal += initialMatrix_[i][neighbor_j] * gaussian1Dkernel_[row];
                }
            }
            tmpMatrix[i].push_back(approxVal);
        }
    }

    // in y direction
    for(unsigned j = 0; j < xn_; ++j) {
        for(unsigned i = 0; i < yn_; ++i) {
            double approxVal = 0.0;
            for(unsigned col = 0; col < kernelSize; ++col) {
                unsigned neighbor_i = i + col - kernelRad_;
                if(neighbor_i >= 0 && neighbor_i < yn_) {
                    approxVal += tmpMatrix[neighbor_i][j] * gaussian1Dkernel_[col];
                }
            }
            initialMatrix_[i][j] = approxVal;
        }
    }
}

I.e., I'm using the same kernel for the boundary elements. I've tested this code on 100x100 matrix and kernel with 2 radius. And, for example, I have quite big difference between elements at 1,97 and 2,97, although in the initial matrix at that positions there is no sharp transition.

Maybe I need to change kernel when calculating approximate values for the boundary elements?

Thanks in advance.

2

2 Answers

1
votes

That’s probably because you are not correctly handling boundary conditions. In your test:

if(neighbor_i >= 0 && neighbor_i < yn_)

the first part will always be true since neighbor_i is an unsigned, hence always positive number. You may want to change that to a signed value slightly modifying its declaration. Your compiler can check for those kind of mistakes for you with the appropriate warning flags (try -Wall -Wextra).

Edit: Actually, the test is probably not the cause of your problem since you are using relatively small images making the value of neighbor_i greater than yn_ when trying to store a negative value inside of it.

Also, please use a library to do the convolutions. In particular there are quite nice and efficient approximations (Canny-Deriche, product in Fourier domain, …) to Gaussian blurs.

0
votes

I solved this issue as below:

Do not normalize the kernel in createGaussianKerenel_1D() function. Instead, do it in a getBlurredThermalMap()function as below:

void getBlurredThermalMap() {
    assert(!gaussian1Dkernel_.empty());
    vector<vector<double> > tmpMatrix(yn_);
    unsigned kernelSize = 2 * kernelRad_ + 1;

    // in x direction
    for(unsigned i = 0; i < yn_; ++i) {
        for(unsigned j = 0; j < xn_; ++j) {
            double approxVal = 0.0;
            double sumNorm = 0.0;
            for(unsigned row = 0; row < kernelSize; ++row) {
                unsigned neighbor_j = j + row - kernelRad_;
                // ignore values that are out of bound
                if(neighbor_j >= 0 && neighbor_j < xn_) {
                    approxVal += initialMatrix_[i][neighbor_j] * gaussian1Dkernel_[row];
                    sumNorm += gaussian1Dkernel_[row];
                }
            }
            approxVal /= sumNorm;
            tmpMatrix[i].push_back(approxVal);
        }
    }

    // in y direction
    for(unsigned j = 0; j < xn_; ++j) {
        for(unsigned i = 0; i < yn_; ++i) {
            double approxVal = 0.0;
            double sumNorm = 0.0;
            for(unsigned col = 0; col < kernelSize; ++col) {
                unsigned neighbor_i = i + col - kernelRad_;
                if(neighbor_i >= 0 && neighbor_i < yn_) {
                    approxVal += tmpMatrix[neighbor_i][j] * gaussian1Dkernel_[col];
                    sumNorm += gaussian1Dkernel_[row];
                }
            }
            approxVal /= sumNorm;
            initialMatrix_[i][j] = approxVal;
        }
    }
}