Last active
December 15, 2022 19:06
-
-
Save nepomnyi/6b95cb0279d9401d8b916ab77a20a3b5 to your computer and use it in GitHub Desktop.
The custom implementation of the PIV image preprocessing algorithm after N. Deen et al., "On image preprocessing for PIV of single- and two-phase flows over reflecting objects", Exp. Fluids, 49(2):525-530, 2010.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import cv2 | |
import numpy as np | |
def customDeen(frame0, | |
frame1, | |
blockOne, | |
algorithmType = 'basic', | |
minMaxKernel = (15,15), | |
minMaxLowLimit = 10, | |
cliplimit = 2.0, | |
tileGridSize = (4, 4) | |
): | |
""" | |
This is the custom implementation of the image processing alogrithm for PIV | |
after Deen, Willems, Annaland, Kuipers, Lammertink, Kemperman, Wessling, | |
Meer, "On image preprocessing for PIV of single and two phase flows over | |
reflecting objects", 2010, Exp Fluids, p.526, Fig.1. | |
I say custom because I substituted the first min-max filtering block with a | |
choise of either adaptive histogram equalization or my custom implementation | |
of min-max filter after Adrian&Westerweel PIV book. | |
Parameters: frame0 (ndarray) - the original image of the first PIV frame | |
frame1 (ndarray) - the original image of the second PIV frame | |
algorithmType (str) - 'basic' which is default or 'advanced'; | |
'basic' is the single phase flow | |
implementation (see Fig.1 in the | |
referenced article), | |
'advanced' is the two phase flow | |
implementation (see Fig.1 in the | |
referenced article) | |
blockOne (str) - what we want to do in block one of Deen's | |
algorithm: | |
either | |
adaptive histogram equalization ("adHist") | |
or | |
normalization using my custom built min-max | |
filter ("minMax") | |
minMaxKernel (ndarray) - the kernel of the min-max filter (see | |
minMaxFilter function) | |
minMaxLowLimit (int) - the low limit imposed on the contrast by | |
the min-max filter (see minMaxFilter | |
function) | |
clipLimit (float) - is the parameter of the adaptive | |
histogram equalization | |
tileGridSize (ndarray) - is the parameter of the adaptive | |
histogram equalization | |
Returns: frame0Subtracted (ndarray) - processed frame0 | |
frame1Subtracted (ndarray) - processed frame1 | |
""" | |
# Right now only 'basic' algorithmType is implemented, thus, no if statements ... | |
# Start with checking if the images are gray scaled | |
if len(frame0.shape) > 2: | |
frame0 = cv2.cvtColor(frame0, cv2.COLOR_BGR2GRAY) | |
if len(frame1.shape) > 2: | |
frame1 = cv2.cvtColor(frame1, cv2.COLOR_BGR2GRAY) | |
# First block: either adaptive histogram equalization or normalization | |
# Adaptive histogram equalization is copied from here: | |
# https://pyimagesearch.com/2021/02/01/opencv-histogram-equalization-and-adaptive-histogram-equalization-clahe/ | |
if blockOne == "adHist": | |
clahe = cv2.createCLAHE(clipLimit = cliplimit, tileGridSize = tileGridSize) | |
frame0Enhanced = clahe.apply(frame0) | |
frame1Enhanced = clahe.apply(frame1) | |
elif blockOne == "minMax": | |
frame0Enhanced = minMaxFilter(frame0, minMaxKernel, minMaxLowLimit) | |
frame1Enhanced = minMaxFilter(frame1, minMaxKernel, minMaxLowLimit) | |
# Second block: stretching | |
frame0Stretched = stretchFilter(frame0Enhanced) | |
frame1Stretched = stretchFilter(frame1Enhanced) | |
# After all these operations, our image is not a 0-255 image anymore. To get | |
# it back to a 0-255 image (an 8-bit integer), do the following (note that | |
# the solution is copied from here: | |
# https://stackoverflow.com/a/61994089/10073233): | |
frame0Stretched = np.uint8(frame0Stretched*255) | |
frame1Stretched = np.uint8(frame1Stretched*255) | |
# Third block: background subtraction. Note, that it can result in negative | |
# values. That's why the formula in the article takes max. | |
# But cv2.subtract takes care of the negative values, so we don't need max. | |
frame0Subtracted = cv2.subtract(frame0Stretched, frame1Stretched) | |
frame1Subtracted = cv2.subtract(frame1Stretched, frame0Stretched) | |
return frame0Subtracted, frame1Subtracted | |
def stretchFilter(img): | |
""" | |
This is the implementation of the stretching step after Deen, Willems, | |
Annaland, Kuipers, Lammertink, Kemperman, Wessling, Meer, "On image | |
preprocessing for PIV of single and two phase flows over reflecting | |
objects", 2010, Exp Fluids, p.526, Fig.1. | |
Parameters: img (ndarray) - original image | |
Returns: stretchedImg (ndarray) - image with the stretched intensities | |
""" | |
# cv2.minMaxLoc finds min and max intensities in an image. It works with | |
# single channel images only, therefore we have to make sure | |
# we don't have more than 1 channel in the original image. I'm not going to | |
# use the if statement for that. Instead, no matter what, automatically | |
# convert the original image into a 1D array using np.flatten(). | |
Imin, Imax, _, _ = cv2.minMaxLoc(img.flatten()) # min and max intensities | |
stretchedImg = (img - Imin) / (Imax - Imin) | |
return stretchedImg | |
def minMaxFilter(img, filterSize, minContrast): | |
""" | |
After R.Adrian, J.Westerweel, "Particle image velocimetry", Cambridge | |
university press, 2011. See ch.6.1.2, p.248-250. | |
Parameters: img (ndarray) - image to be filtered | |
filterSize (ndarray) - a 1x2 numpy array of the filter height | |
and width correspondingly | |
minContrast (float) - minimum contrast value imposed on the | |
image (if the calculated contrast falls | |
bellow this level, this level is imposed | |
as the minimum contrast) | |
Returns: imgFiltered (ndarray) - filtered image | |
""" | |
# Define the lower and upper envelopes (see the book) of the image intensity | |
low = cv2.erode(img, | |
cv2.getStructuringElement(cv2.MORPH_RECT, filterSize)) | |
upp = cv2.dilate(img, | |
cv2.getStructuringElement(cv2.MORPH_RECT, filterSize)) | |
# Smooth the lower and upper envelopes with a uniform filter of the same size | |
low = lowPassFilter(low, filterSize) | |
upp = lowPassFilter(upp, filterSize) | |
# Define contrast and put a lower limit on it | |
contrast = cv2.subtract(upp, low) | |
contrast[contrast < minContrast] = minContrast | |
# Normalize image intensity | |
imgFiltered = np.uint8(cv2.divide(cv2.subtract(img, low), contrast) * 255) | |
return imgFiltered | |
def lowPassFilter(img_src, kernelSize): | |
""" | |
This is a uniform low pass filter, which measn that all the values in the | |
filter kernel are equal to 1. | |
Parameters: img_src (ndarray) - image to be filtered | |
kernelSize (ndarray) - 1x2 numpy array where the first value is | |
the kernel's height, the second value is | |
the kernel's width | |
Returns: img_rst (ndarray) - filtered image | |
""" | |
kernel = np.ones(kernelSize) | |
kernel = kernel/(np.sum(kernel) if np.sum(kernel)!=0 else 1) | |
img_rst = cv2.filter2D(img_src, -1, kernel) | |
return img_rst |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
For the first block of Deen's algorithm I give two options to choose from: either adaptive histogram equalization or min-max filtering after Adrian&Westerweel. That's because sometimes histogram equalization performs better and pretty much always it is easier to start with (with min-max filter one has to tune the parameters before one gets a nice image, whereas histogram equalization gives you a nice image at once and then you have to tune the parameters just slightly). The both methods fall within the category of image normalization.
With respect to the resultant velocity fields, the both methods yield equal results.
Also, note that I didn't implement min-max filter as given in Deen's article. Instead, I implemented it as given in Adrian&Westerweel. Here's my separate Gist with the min-max filter implementation.