Mejdi Dallel
Mejdi Dallel

Reputation: 642

How to find the Gaussian Weighted-Average and Standard deviation of a structural element

I'm trying to implement an Intensity Normalization algorithm that is described by this formula:

x' = (x - gaussian_weighted_average) / std_deviation

The paper I'm following describes that I have to find the gaussian weighted average and the standard deviation corresponding to each pixel "x" neighbors using a 7x7 kernel.

PS: x' is the new pixel value.

enter image description here

So, my question is: how can I compute a gaussian weighted average and the standard deviation for each pixel in image using a 7x7 kernel?

Does OpenCV provide any method to solve this?

import cv2
img = cv2.imread("b.png", 0)
widht = img.shape[0]
height = img.shape[1]
for i in range (widht):
    for j in range (height):
        new_image = np.zeros((height,width,1), np.uint8)
        new_image[i][j] = img[i][j] - ...

Upvotes: 1

Views: 2195

Answers (1)

Berriel
Berriel

Reputation: 13641

The original implementation (C++) of the author can be found here: see GenerateIntensityNormalizedDatabase().

This has been re-implemented by another student in python. The python implementation is:

import cv2
import numpy as np

def StdDev(img, meanPoint, point, kSize):
    kSizeX, kSizeY = kSize / 2, kSize / 2

    ystart = point[1] - kSizeY if 0 < point[1] - kSizeY < img.shape[0] else 0
    yend = point[1] + kSizeY + 1 if 0 < point[1] + kSizeY + 1 < img.shape[0] else img.shape[0] - 1

    xstart = point[0] - kSizeX if 0 < point[0] - kSizeX < img.shape[1] else 0
    xend = point[0] + kSizeX + 1 if 0 < point[0] + kSizeX + 1 < img.shape[1] else img.shape[1] - 1

    patch = (img[ystart:yend, xstart:xend] - meanPoint) ** 2
    total = np.sum(patch)
    n = patch.size

    return 1 if total == 0 or n == 0 else np.sqrt(total / float(n))


def IntensityNormalization(img, kSize):
    blur = cv2.GaussianBlur(img, (kSize, kSize), 0, 0).astype(np.float64)
    newImg = np.ones(img.shape, dtype=np.float64) * 127

    for x in range(img.shape[1]):
        for y in range(img.shape[0]):
            original = img[y, x]
            gauss = blur[y, x]
            desvio = StdDev(img, gauss, [x, y], kSize)

            novoPixel = 127
            if desvio > 0:
                novoPixel = (original - gauss) / float(desvio)

            newVal = np.clip((novoPixel * 127 / float(2.0)) + 127, 0, 255)
            newImg[y, x] = newVal
    return newImg

To use the intensity normalization, you could do this:

kSize = 7
img = cv2.imread('{IMG_FILENAME}', cv2.IMREAD_GRAYSCALE).astype(np.float64)
out = IntensityNormalization(img, kSize)

To visualize the resulting image, don't forget to convert out back to np.uint8 (why?). I'd recommend you to use the original implementation in C++ if you want to reproduce his results.

Disclaimer: I'm from the same lab of the author of this paper.

Upvotes: 2

Related Questions