Issues in Gaussian blur transform with kernel computation and normalization

I am currently trying to implement Gaussian blurring myself for learning purposes. The blurring effect does not take place. Instead the image contrast is enhanced.

#include <algorithm>
#include <math.h>
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/opencv.hpp>
#include <vector>
using namespace std;
using namespace cv;

void displayImage(Mat &img, unsigned int time = 0, string title = "frame") {
  imshow(title, img);
  waitKey(time);
}

bool isValidPoint(Mat &img, int x, int y) {
  int rows = img.rows;
  int cols = img.cols;
  return (x >= 0 and x < rows and y >= 0 and y < cols);
}

static double square(double x) { return x * x; }

static double computeGaussianFunc(double x, double y, double mu, double sigma) {
  double val =
      exp(-0.5 * ((square((x - mu) / sigma)) + square((y - mu) / sigma))) /
      (2 * M_PI * sigma * sigma);
  return val;
}

double norm2(double val,double minVal,double maxVal){
  double range = (maxVal-minVal);
  double newVal  = (val-minVal)/range;
  return (int(255*newVal));
}

//  oldMin,  oldMax,  newMin,  newMax,  oldVal
double getNormalizedValue(double oldMin, double oldMax, double newMin,
                          double newMax, double oldVal) {
  double oldRange = (oldMax - oldMin);
  double newRange = (newMax - newMin);
  double newVal = (newMin + ((newRange * (oldVal - oldMin)) / (oldRange)));
  return newVal;
}

static vector<vector<double>> getGuassianKernal(double sigma) {
  int size = 2 * ceil(3 * sigma) + 1;
  vector<vector<double>> Kernal(size, vector<double>(size, 0.0));
  double sum = 0.0;
  int center = size / 2;
  for (int i = 0; i < size; i++) {
    for (int j = 0; j < size; j++) {
      int x = i - center;
      int y = j - center;
      Kernal[i][j] = computeGaussianFunc(i, j, size / 2, sigma);
      sum += Kernal[i][j];
    }
  }
  if (sum != 0) {
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        Kernal[i][j] /= sum;
      }
    }
  }
  return Kernal;
}

static Mat guassianFilterTransformExp(Mat &img, double Sigma) {
  vector<vector<double>> filter = getGuassianKernal(Sigma);
  int FiltSize = filter.size();
  int trows = img.rows - FiltSize + 1;
  int tcols = img.cols - FiltSize + 1;
  // Final output
  Mat transformed(trows, tcols, CV_8U);
  // intermidiate matrix using for normalizing.
  Mat inter(trows, tcols, CV_64F);
  // min and max values of value after convolving with filter for normalization
  double minVal = 10000.0;
  double maxVal = 0.0;
  for (int i = 1; i <= trows - 2; i++) {
    for (int j = 1; j <= tcols - 2; j++) {
      double tval = 0;
      for (int x = -1; x <= 1; x++) {
        for (int y = -1; y <= 1; y++) {
          tval = tval + (filter[x + 1][y + 1] *
                         static_cast<double>(img.at<u_char>(i + x, j + y)));
          minVal = min(minVal,tval);
          maxVal = max(maxVal,tval);
        }
      }
      inter.at<double>(i,j) = tval;
    }
  }
  for (int i = 0; i < trows; i++) {
    for (int j = 0; j <tcols; j++) {
      double val = inter.at<double>(i,j);
      double newVal = norm2(val,minVal,maxVal);
      cout<<"Computed Val : "<<val<<" Src : "<<static_cast<double>(img.at<u_char>(i, j))<<"  Normalized : "<<newVal<<'\n';
      transformed.at<u_char>(i,j) = static_cast<u_char>(min(255.0,newVal));
    }
  }
  return (transformed);
}


int main() {
  string imPath =
      "/home/panirpal/workspace/Projects/ComputerVision/data/images/chess2.jpg";
  Mat img = imread(imPath, IMREAD_GRAYSCALE);
  if (!img.empty()) {
    displayImage(img);
    Mat out = guassianFilterTransformExp(img,1.3);
    displayImage(out);
  } else
    cerr << "image not found! exiting...";
  return 0;
}

On debugging I get multiple observations. First of all, I am computing kernel size based on sigma 2 * ceil(3 * sigma) + 1; found on the Internet. So kernel size comes to be 9 ( sigma given = 1.3 , 2*(ceil(3*1.3) + 1 ). And kernel values come to be pretty small, when the filter is convolved/cross co-relation operation with the image patch, the resulting value ( denoted as tval in code) comes to be fairly small, and I suppose no surprises there. Then I try to normalize this tval to a range 0-255, since rendering tval as is would be resulting in a near blacked out image, I used one of the normalization techniques mentioned here and following is some logs I collected to get a sense of what is going on, Computed val stands for tval Src stands for the orignal ( grayscale pixel value in source image ) and Normalized stands for newVal in code. The strange part is the Normalized values are higher than source they should be lower and smoothened. And the output also reflects the same. What could I be doing wrong here?

Computed Val : 2.25082 Src : 167  Normalized : 177
Computed Val : 2.21529 Src : 161  Normalized : 174
Computed Val : 2.27372 Src : 159  Normalized : 179
Computed Val : 2.36755 Src : 160  Normalized : 186
Computed Val : 2.39062 Src : 177  Normalized : 188
Computed Val : 2.32093 Src : 174  Normalized : 182
Computed Val : 1.25675 Src : 168  Normalized : 99
Computed Val : 0.495057 Src : 91  Normalized : 39
Computed Val : 1.14986 Src : 39  Normalized : 90
Computed Val : 2.07179 Src : 64  Normalized : 163
Computed Val : 2.34365 Src : 132  Normalized : 184
Computed Val : 2.39052 Src : 172  Normalized : 188

Input image (named chess2.jpg):

Output:

Note: ensure OpenCV is built on your machine. Can be found in the OpenCV install guide.

The command used to build the code:

g++ exp.cpp -I/usr/local/include/opencv4 -Wl,-rpath,/usr/local/lib /usr/local/lib/libopencv_highgui.so.4.8.0 /usr/local/lib/libopencv_ml.so.4.8.0 /usr/local/lib/libopencv_objdetect.so.4.8.0 /usr/local/lib/libopencv_photo.so.4.8.0 /usr/local/lib/libopencv_stitching.so.4.8.0 /usr/local/lib/libopencv_video.so.4.8.0 /usr/local/lib/libopencv_videoio.so.4.8.0 /usr/local/lib/libopencv_imgcodecs.so.4.8.0 /usr/local/lib/libopencv_calib3d.so.4.8.0 /usr/local/lib/libopencv_dnn.so.4.8.0 /usr/local/lib/libopencv_features2d.so.4.8.0 /usr/local/lib/libopencv_flann.so.4.8.0 /usr/local/lib/libopencv_imgproc.so.4.8.0 /usr/local/lib/libopencv_core.so.4.8.0 -lm -o exp

Expected output to be blurred image.

In your Gaussian filter function, you create a Gaussian kernel of the right sizes, you read those sizes and adjust the loop over the image according to it, but then the inner loops are:

      for (int x = -1; x <= 1; x++) {
        for (int y = -1; y <= 1; y++) {

which only iterate over the top-left 3x3 region of the kernel.

I don't like it when a filter function returns a smaller image than it got as input, because it makes it hard to combine the filter result with the original. But given that you do return a smaller image, you could simplify your code by not worrying about where the center of the kernel is:

static Mat guassianFilterTransformExp(Mat &img, double Sigma) {
  vector<vector<double>> filter = getGuassianKernal(Sigma);
  int FiltSize = filter.size();
  int trows = img.rows - FiltSize + 1;
  int tcols = img.cols - FiltSize + 1;
  Mat transformed(trows, tcols, CV_8U);
  for (int i = 0; i < trows; i++) {
    for (int j = 0; j < tcols; j++) {
      double tval = 0;
      for (int x = 0; x < FiltSize; x++) {
        for (int y = 0; y < FiltSize; y++) {
          tval = tval + (filter[x][y] *
                         static_cast<double>(img.at<u_char>(i + x, j + y)));
        }
      }
      tval = std::min(tval, 255.0);
      tval = std::max(tval, 0.0);
      transformed.at<u_char>(i,j) = static_cast<u_char>(std::round(tval));
    }
  }
  return (transformed);
}

This simply removes all the +1 and -1 operations, which were meant for the 3x3 kernel and would have been replaced by +FiltSize/2 and -FiltSize/2. (I've also removed the normalization part, which you shouldn't need.)

Note that the Gaussian is separable, which means you can implement this filter much more efficiently. Also vector<vector<double>> is a highly inefficient way to store a 2D matrix, and img.at<u_char>() is the least efficient way to access a pixel when in a loop over all image pixels. But I assume this is a learning exercise and not meant for serious production work. If this is meant for production, please use a library that already implements the operation efficiently.

There are other issues with the code, I didn't comb though all of it. For example, as Marek R mentioned in a comment,

Kernal[i][j] = computeGaussianFunc(i, j, size / 2, sigma);

should be

Kernal[i][j] = computeGaussianFunc(x, y, size / 2, sigma);
                                   ^^^^^