Image goes black after thresholding

Question

I am trying to extract blood network from this face image: Face image

For such task, i am using the P&M anisotropic diffusion found in this question: Anisotropic diffusion 2d images. Then i am using tophat transform followed by blackhat transform, afterwards i use a simple threshold to set to 255 all pixel that has an intensity value of 100.

The problem is that, after i use the threshold and try to open the image, whatever way i try, the image is displayed as fully black:

In short, my goal is to extract the blood vessels using P&M anisotropic diffusion with structuring element of flat disk of 5x5, then apply tophat and blackhat, respectively and a simple threshold and actually be able to view the image afterwards.

Here's my code on how i am trying it:

    import cv2 
    import import cv2 numpy as np
    import warnings

    face_img=mpimg.imread('path')
    def anisodiff(img, niter=1, kappa=50, gamma=0.1, step=(1., 1.), option=1):
        if img.ndim == 3:
            m = "Only grayscale images allowed, converting to 2D matrix"
            warnings.warn(m)
            img = img.mean(2)
    
        img = img.astype('float32')
        imgout = img.copy()
    
        deltaS = np.zeros_like(imgout)
        deltaE = deltaS.copy()
        NS = deltaS.copy()
        EW = deltaS.copy()
        gS = np.ones_like(imgout)
        gE = gS.copy()
    
        for ii in range(niter):
    
            deltaS[:-1, :] = np.diff(imgout, axis=0)
            deltaE[:, :-1] = np.diff(imgout, axis=1)
    
            if option == 1:
                    gS = np.exp(-(deltaS/kappa)**2.)/step[0]
                    gE = np.exp(-(deltaE/kappa)**2.)/step[1]
            elif option == 2:
                    gS = 1./(1.+(deltaS/kappa)**2.)/step[0]
                    gE = 1./(1.+(deltaE/kappa)**2.)/step[1]
    
            E = gE*deltaE
            S = gS*deltaS
    
            NS[:] = S
            EW[:] = E
            NS[1:, :] -= S[:-1, :]
            EW[:, 1:] -= E[:, :-1]
    
            imgout += gamma*(NS+EW)
    
        return imgout

new_img = anisodiff(face_img, niter=1, kappa=20, gamma=0.1, step=(1., 1.), option=1)

filterSize =(3, 3) 
kernel = cv2.getStructuringElement(cv2.MORPH_RECT,  
                                   filterSize) 
input_image = new_img 
first_tophat_img = cv2.morphologyEx(input_image,  
                              cv2.MORPH_TOPHAT, 
                              kernel)
filterSize =(3, 3) 
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, 
                                   filterSize) 
  
second_tophat_img = cv2.morphologyEx(input_image,  
                              cv2.MORPH_BLACKHAT, 
                              kernel) 

ret, thresh1 = cv2.threshold(second_tophat_img, 200, 255, cv2.THRESH_BINARY)

Even when i set the threshold to 254 for instance, the image goes black.

user_4926757 · Answer

I executed a simple MATLAB implementation, and got a nice result.

MATLAB code:

I = imread('02_giorgos_1_f_M_30_830.tif');
I = im2double(uint8(I));

J = imdiffusefilt(I);
K = imtophat(J, ones(3));

figure;imshow(imadjust(K, stretchlim(K)));

Result:

I don't know if you know MATLAB, but I used the default arguments of imdiffusefilt (equivalent to anisodiff in your code).

Default MATLAB arguments are equivalent to:

Input image is in range [0, 1] and not [0, 255].
niter=5 (note: you used only 1 iteration and it's not enough).
kappa=0.1
gamma=0.125
MATLAB default is 8 neighbors connectivity (not 4 neighbors as used in anisodiff).

8 neighbors connectivity:

For getting same result as in MATLAB, I implemented an 8 neighbors connectivity Anisotropic diffusion (based on MATLAB source code).
Note: with 4 neighbors connectivity it's working, but result is not so nice as using 8 neighbors.

Displaying the output image:

In order to display the output image correctly, I used imadjust(K, stretchlim(K)).
The command stretches the range of the input image such that percentile 1 goes to 0, and percentile 99 goes to 1 (linear stretch).

One more thing:

Instead of using fixed threshold of 200, I used percentile 95 threshold:

 t = np.percentile(first_tophat_img, 95)
 ret, thresh1 = cv2.threshold(first_tophat_img, t, 255, 
 cv2.THRESH_BINARY)

Here is the code (uses cv2.imshow for testing):

import cv2
import numpy as np
import matplotlib.image as mpimg
import warnings

face_img = mpimg.imread('02_giorgos_1_f_M_30_830.tif')

def anisodiff8neighbors(img, niter=5, kappa=0.1, gamma=0.125):
    """ See https://www.mathworks.com/help/images/ref/imdiffusefilt.html 
        Anisotropic diffusion filtering with 8 neighbors
        Range of img is assumed to be [0, 1] (not [0, 255]).
    """
    if img.ndim == 3:
        m = "Only grayscale images allowed, converting to 2D matrix"
        warnings.warn(m)
        img = img.mean(2)
    
    img = img.astype('float32')
    imgout = img.copy()
    
   
    for ii in range(niter):
        # MATLAB source code is commented
        #paddedImg = padarray(I, [1 1], 'replicate');
        padded_img = np.pad(imgout, (1, 1), 'edge')

        #diffImgNorth = paddedImg(1:end-1,2:end-1) - paddedImg(2:end,2:end-1);
        #diffImgEast = paddedImg(2:end-1,2:end) - paddedImg(2:end-1,1:end-1);
        #diffImgNorthWest = paddedImg(1:end-2,1:end-2) - I;
        #diffImgNorthEast = paddedImg(1:end-2,3:end) - I;
        #diffImgSouthWest = paddedImg(3:end,1:end-2) - I;
        #diffImgSouthEast = paddedImg(3:end,3:end) - I;
        diff_img_north = padded_img[0:-1, 1:-1] - padded_img[1:, 1:-1]
        diff_img_east = padded_img[1:-1, 1:] - padded_img[1:-1, 0:-1]
        diff_img_north_west = padded_img[0:-2, 0:-2] - imgout
        diff_img_north_east = padded_img[0:-2, 2:] - imgout
        diff_img_south_west = padded_img[2:, 0:-2] - imgout
        diff_img_south_east = padded_img[2:, 2:] - imgout

        #case 'exponential'
        #conductCoeffNorth = exp(-(abs(diffImgNorth)/gradientThreshold).^2);
        #conductCoeffEast = exp(-(abs(diffImgEast)/gradientThreshold).^2);
        #conductCoeffNorthWest = exp(-(abs(diffImgNorthWest)/gradientThreshold).^2);
        #conductCoeffNorthEast = exp(-(abs(diffImgNorthEast)/gradientThreshold).^2);
        #conductCoeffSouthWest = exp(-(abs(diffImgSouthWest)/gradientThreshold).^2);
        #conductCoeffSouthEast = exp(-(abs(diffImgSouthEast)/gradientThreshold).^2);
        conduct_coeff_north = np.exp(-(np.abs(diff_img_north)/kappa)**2.0)
        conduct_coeff_east = np.exp(-(np.abs(diff_img_east)/kappa)**2.0)
        conduct_coeff_north_west = np.exp(-(np.abs(diff_img_north_west)/kappa)**2.0)
        conduct_coeff_north_east = np.exp(-(np.abs(diff_img_north_east)/kappa)**2.0)
        conduct_coeff_south_west = np.exp(-(np.abs(diff_img_south_west)/kappa)**2.0)
        conduct_coeff_south_east = np.exp(-(np.abs(diff_img_south_east)/kappa)**2.0)
        
        #fluxNorth = conductCoeffNorth .* diffImgNorth;
        #fluxEast =  conductCoeffEast .* diffImgEast;
        #fluxNorthWest = conductCoeffNorthWest .* diffImgNorthWest;
        #fluxNorthEast = conductCoeffNorthEast .* diffImgNorthEast;
        #fluxSouthWest = conductCoeffSouthWest .* diffImgSouthWest;
        #fluxSouthEast = conductCoeffSouthEast .* diffImgSouthEast;
        flux_north = conduct_coeff_north * diff_img_north
        flux_east = conduct_coeff_east * diff_img_east
        flux_north_west = conduct_coeff_north_west * diff_img_north_west
        flux_north_east = conduct_coeff_north_east * diff_img_north_east
        flux_south_west = conduct_coeff_south_west * diff_img_south_west
        flux_south_east = conduct_coeff_south_east * diff_img_south_east

        #% Discrete PDE solution
        #I = I + diffusionRate * (fluxNorth(1:end-1,:) - fluxNorth(2:end,:) + ...
        #    fluxEast(:,2:end) - fluxEast(:,1:end-1) + (1/(dd^2)).* fluxNorthWest + ...
        #    (1/(dd^2)).* fluxNorthEast + (1/(dd^2)).* fluxSouthWest + (1/(dd^2)).* fluxSouthEast);
        imgout = imgout + gamma * (flux_north[0:-1,:] - flux_north[1:,:] + 
            flux_east[:,1:] - flux_east[:,0:-1] + 0.5*flux_north_west + 
            0.5*flux_north_east + 0.5*flux_south_west + 0.5*flux_south_east)
    
    return imgout


#new_img = anisodiff(face_img, niter=1, kappa=20, gamma=0.1, step=(1., 1.), option=1)
face_img = face_img.astype(float) / 255;
#new_img = anisodiff(face_img, niter=5, kappa=0.1, gamma=0.125, step=(1., 1.), option=1)
new_img = anisodiff8neighbors(face_img, niter=5, kappa=0.1, gamma=0.125)

filterSize =(3, 3) 
kernel = cv2.getStructuringElement(cv2.MORPH_RECT,  
                                   filterSize) 
input_image = new_img 
first_tophat_img = cv2.morphologyEx(input_image,
                              cv2.MORPH_TOPHAT, 
                              kernel)

# Use percentile 95 (of image) as threshold instead of fixed threshold 200
t = np.percentile(first_tophat_img, 95)
ret, thresh1 = cv2.threshold(first_tophat_img, t, 255, cv2.THRESH_BINARY)
cv2.imshow('thresh1', thresh1)


filterSize =(3, 3) 
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, 
                                   filterSize) 
  

second_tophat_img = cv2.morphologyEx(input_image,  
                              cv2.MORPH_BLACKHAT, 
                              kernel) 


#ret, thresh1 = cv2.threshold(second_tophat_img, 200, 255, cv2.THRESH_BINARY)

# Use percentile 95 (of image) as threshold instead of fixed threshold 200
t = np.percentile(second_tophat_img, 95)
ret, thresh2 = cv2.threshold(second_tophat_img, t, 255, cv2.THRESH_BINARY)
cv2.imshow('thresh2', thresh2)

lo, hi = np.percentile(first_tophat_img, (1, 99))
first_tophat_img_stretched = (first_tophat_img.astype(float) - lo) / (hi-lo) # Apply linear "stretch" - lo goes to 0, and hi goes to 1
cv2.imshow('first_tophat_img_stretched', first_tophat_img_stretched)


cv2.waitKey()
cv2.destroyAllWindows()

Result:

thresh1:

thresh2:

first_tophat_img_stretched:

Image goes black after thresholding

1 Answers