imflatfield MATLAB for Python use

Question

I am trying to find an equivalent Python function for MATLAB imflatfield function.

I have a section of code that modifies an image and I want to convert it to Python.

Here is the MATLAB code:

I = imread('lcs2.png');

out2 = imflatfield(I,30);

shadow_lab = rgb2lab(out2);
max_luminosity = 100;
L = shadow_lab(:,:,1)/max_luminosity;
shadow_adapthisteq = shadow_lab;
shadow_adapthisteq(:,:,1) = adapthisteq(L)*max_luminosity;
shadow_adapthisteq = lab2rgb(shadow_adapthisteq);

imwrite(shadow_adapthisteq,'lcs2_adap.jpg');

Original image

Final results from MATLAB

Rotem · Answer 1 · 2020-04-08T20:18:30.137

Since MATLAB releases the source code of imflatfield, it is not so difficult to implement it in Python using OpenCV.

Note: The implementation is specific to uint8 type and colored image (BGR format in Python).

Here is a MATLAB "manual" implementation of imflatfield:

function B = my_imflatfield(I, sigma)
    A = im2single(I);
    Ihsv = rgb2hsv(A);
    A = Ihsv(:,:,3);

    filterSize = 2*ceil(2*sigma)+1;

    shading = imgaussfilt(A, sigma, 'Padding', 'symmetric', 'FilterSize', filterSize); % Calculate shading
    meanVal = mean(A(:),'omitnan');

    % Limit minimum to 1e-6 instead of testing using isnan and isinf after division.
    shading = max(shading, 1e-6);

    B = A*meanVal./shading;

    %B(isnan(B)) = 0; % sometimes instances of 0/0 happen, making NaN values.
    %B(isinf(B)) = 0; % sometimes values are divided by 0, making Inf values.

    % Put processed V channel back into HSV image, convert to RGB
    Ihsv(:,:,3) = B;
    B = hsv2rgb(Ihsv);

    B = im2uint8(B);
end

Here is an equivalent Python implementation (using OpenCV):

import cv2
import numpy as np

def imflatfield(I, sigma):
    """Python equivalent imflatfield implementation
       I format must be BGR and type of I must be uint8"""
    A = I.astype(np.float32) / 255  # A = im2single(I);
    Ihsv = cv2.cvtColor(A, cv2.COLOR_BGR2HSV)  # Ihsv = rgb2hsv(A);
    A = Ihsv[:, :, 2]  # A = Ihsv(:,:,3);

    filterSize = int(2*np.ceil(2*sigma) + 1);  # filterSize = 2*ceil(2*sigma)+1;

    # shading = imgaussfilt(A, sigma, 'Padding', 'symmetric', 'FilterSize', filterSize); % Calculate shading
    shading = cv2.GaussianBlur(A, (filterSize, filterSize), sigma, borderType=cv2.BORDER_REFLECT)

    meanVal = np.mean(A)  # meanVal = mean(A(:),'omitnan')

    #% Limit minimum to 1e-6 instead of testing using isnan and isinf after division.
    shading = np.maximum(shading, 1e-6)  # shading = max(shading, 1e-6);

    B = A*meanVal / shading  # B = A*meanVal./shading;

    #% Put processed V channel back into HSV image, convert to RGB
    Ihsv[:, :, 2] = B  # Ihsv(:,:,3) = B;

    B = cv2.cvtColor(Ihsv, cv2.COLOR_HSV2BGR)  # B = hsv2rgb(Ihsv);

    B = np.round(np.clip(B*255, 0, 255)).astype(np.uint8)  # B = im2uint8(B);

    return B

# Read input imgae
I = cv2.imread('destroyer.jpg')

sigma = 30

out2 = imflatfield(I, sigma)

cv2.imwrite('imflatfield_py_destroyer.png', out2)

The above implementation reads the input image, and write the result to image file.

Comparing results using MATLAB (for testing):

I = imread('destroyer.jpg');

out1 = imflatfield(I, 30);

out2 = my_imflatfield(I, 30);

% Compare results of imflatfield and my_imflatfield:
max(max(max(imabsdiff(out1, out2))))

% figure;imshow(out2)

imwrite(out2, 'imflatfield_destroyer.png');

% Read Python result
out3 = imread('imflatfield_py_destroyer.png');

% Compare results of imflatfield and Python imflatfield:
max(max(max(imabsdiff(out1, out3))))

The maximum absolute difference between MATALB imflatfield and my_imflatfield is 0.
The maximum absolute difference between MATALB imflatfield and Python imflatfield is 1.

Converting the complete MATLAB code to Python:

sigma = 30

out2 = imflatfield(I, sigma)


# Conver out2 to float32 before converting to LAB
out2 = out2.astype(np.float32) / 255  # out2 = im2single(out2);

shadow_lab = cv2.cvtColor(out2, cv2.COLOR_BGR2Lab)  # shadow_lab = rgb2lab(out2);

max_luminosity = 100
L = shadow_lab[:, :, 0] / max_luminosity  # L = shadow_lab(:,:,1)/max_luminosity;


shadow_adapthisteq = shadow_lab.copy() # shadow_adapthisteq = shadow_lab;

# shadow_adapthisteq(:,:,1) = adapthisteq(L)*max_luminosity;
clahe = cv2.createCLAHE(clipLimit=20, tileGridSize=(8,8))
cl1 = clahe.apply((L*(2**16-1)).astype(np.uint16)) # CLAHE in OpenCV does not support float32 (convert to uint16 and back).
shadow_adapthisteq[:, :, 0] = cl1.astype(np.float32) * max_luminosity / (2**16-1)

shadow_adapthisteq = cv2.cvtColor(shadow_adapthisteq, cv2.COLOR_Lab2BGR) # shadow_adapthisteq = lab2rgb(shadow_adapthisteq);

# Convert shadow_adapthisteq to uint8
shadow_adapthisteq = np.round(np.clip(shadow_adapthisteq*255, 0, 255)).astype(np.uint8)  # B = im2uint8(B);

cv2.imwrite('shadow_adapthisteq.jpg', shadow_adapthisteq)  # imwrite(shadow_adapthisteq,'lcs2_adap.jpg');

Result is not going to be identical to MATLAB, because adapthisteq in MATLAB is not identical to CLAHE in OpenCV.

Result:

imflatfield MATLAB for Python use

1 Answers1

Linked