Correct way to compute coefficient from wavelet decomposition (Need Explanation)

[ghost]

Hello, im trying to do mammogram image enhancement in wavelet domain based on this paper A Direct Image Contrast Enhancement Algorithm in the Wavelet Domain for Screening Mammograms where its used 4 level decomposition with 'db16' wavelet. The proposed algorithm works by doing some computation on each level coefficient. For easier and faster reading here i provided a screenshoot of the paper describing the algorithm : Following the algorithm here is my code (based on my understanding of the algorithm itself) :

LAMBDA = 1.80

#based on Eq.10 in the image
#coefficient computation
def compute(L,a,b,c):
    a_b = np.divide(a,b,out=np.zeros_like(a),where = b!=0)
    result = L * a_b * c
    print('max',np.max(result))
    return result

# For 4-level decomposition, 
C = pywt.wavedec2(img, 'db16', mode='periodization', level=4)
# C =  [cA4,(cH4,cV4,cD4), (cH3,cV3,cD3),(cH2,cV2,cD2),(cH1,cV1,cD1)]
# C =  (C[0],   (C[-4]),      (C[-3]),      (C[-2]),      (C[-1]))

#--------Eq8 here

cA4 = C[0]
lcA4 = cA4
# for cA3,
C3=(cA4,(C[-4]))
cA3=pywt.waverec2(C3, 'db16', mode='periodization')
#the computation
print(1)
lcD4 = tuple([LAMBDA * v for v in C[-4]])
LC3 = (lcA4,lcD4)
lcA3 = pywt.waverec2(LC3,'db16',mode='periodization')

#--------Eq10 here

# for cA2,
C2 = (cA3,(C[-3]))
cA2 = pywt.waverec2(C2,'db16',mode='periodization')
#the computation
print(2)
lcD3 = tuple([compute(LAMBDA,v,cA3,lcA3) for v in C[-3]])
LC2 = (lcA3,lcD3)
lcA2 = pywt.waverec2(LC2,'db16',mode='periodization')

#for cA1,
C1 = (cA2,(C[-2]))
cA1 = pywt.waverec2(C1,'db16',mode='periodization')
#the computation
print(3)
lcD2 = tuple([compute(LAMBDA,v,cA2,lcA2) for v in C[-2]])
LC1 = (lcA2,lcD2)
lcA1 = pywt.waverec2(LC1,'db16',mode='periodization')

#for cA0
C0 = (cA1,(C[-1]))
cA0 = pywt.waverec2(C0,'db16',mode='periodization')
# the computation
print(4)
lcD1 = tuple([compute(LAMBDA,v,cA1,lcA1) for v in C[-1]])
LC0 = (lcA1,lcD1)
lcA0 = pywt.waverec2(LC0,'db16',mode='periodization')

enhanced_img = lcA0.astype(np.uint8)

This is the original image (a grayscale mammogram image) that im trying to enhance: The original image have resolution of 1024 by 1024, i keep the image size 2^n by 2^n to avoid size error while reconstructing back. I used lettebox resizing method (to keep the aspect ratio of the image) before wavelet decomposing. This image has maximum pixel value of 242 and minimum pixel value of 0.

The result of the enhanced image (lcA0) is this: As you can see there is alot of distortion in area that supposed to be black while (i think) there is enhancement in breast region. FYI the image before casted to uint8 has maximum pixel value of 9.139444144535558e+18 and minimum pixel value of -6.3183869251787e+18 and the view is all white in breast region with same distortion. After casted to uint8 the maximum pixel value is 255 with minimum of 0.

My hypothesis for the cause of this is due to the computation that i did in those wavelet coefficient as the max value of the resulted computation is so big like 1e+16 and more. But I dont know if its correct or not to point this as the cause. I also handled zero divison that might occur in the computation by setting the result to 0 and im not sure too if its the right way to do it.

Im asking because there are no many example of wavelet coefficient manipulation/computation and the value of the coefficient itself is kinda arbitary for me. So i need some explanation to handle the value of computed coefficient if i ever needed or the right way to recostruct the image using the computed coefficient.Or maybe there other problem causing this.

Sorry this is rather a very long post but im not issuing any problem i just need your help to implement this enhancement techniques. Thankyou.

[grlee77]

Feel free to continue the discussion here, but I am going to close the issue as this page is primarily to report bugs or request specific features for the software rather than for general technical support.

We do have a mailing list for general discussion / technical questions, which while not super active does occasionally get community responses to questions.

[grlee77]

You might want to try normalizing the image roughly to the range [-1, 1] first with something like:

image = image / np.abs(image).max()

I wouldn't cast to uint8 as PyWavelets is going to use floating point anyway for the transforms.

This looks like you need to define a foreground mask (i.e. True where there is signal and zero for background regions). One potential tool for auto-selecting a threshold is Otsu thresholding, but you may have luck just smoothing the image and then adjusting a threshold manually until you get a good mask.

One you have the mask you can then set all values not in the mask to zero after your coefficient division operations:

image[~mask] = 0

You can also potentially avoid divide by zero in the first place by just adding a very small constant to the denominator (i.e. add 1e-6 to the denominator values).

[ghost]

Hi @grlee77 thanks for the kind answer. Im sorry for opening this as an issue, never meant to, later on i will use the mailing list. Before decomposing i cast the image to float32 then only after reconstructing the final coefficient (cA0) i cast it back to uint8. For normalizing the image, i already tried this but i got different result, what i mean is when recostructing the final coefficient then i normalized it back to the initial range after that cast it to uint8 but the image became so dark. Even when im not normalized it back still got dark image. By using a mask what you mean is to mask out all the distorted areas ? For division by zero that actually a nice suggestion i will try it. Thank you so much sir.