Reputation: 11
I want to apply a Gaussian filter on an float numpy array using python2.
In fact I'm trying to rewrite the code Retrospective Correction using Homomorphic Filtering in python,
g(x,y) = exp(LPF(log(f(x,y)))) . C
where LPF(f(x,y)) is the low-pass filter of f(x,y) and C is the normalization coefficient.
Here is my code, which would crash at line:LPFF = scipy_filters.gaussian_filter1d(logF, sigma=10)
fxy = greyscale_matrix
logF = np.log1p(fxy)
print logF
LPFF = scipy_filters.gaussian_filter1d(logF, sigma=10)
print LPFF
expF = np.expm1(LPFF)
tmp_div = fxy / expF
meanF_down = np.mean(tmp_div)
meanF_top = np.mean(fxy)
C = meanF_top / meanF_down
gx_matrix = expF * C
rtn_matrix = greyscale_matrix - gx_matrix
Image.fromarray(gx_matrix).show()
Image.fromarray(rtn_matrix).show()
I try to debug it, when going to
def correlate1d(*args, **kwargs): # real signature unknown
pass
from _nd_image.py the code seems crash with no idea. I also find a github doing this job as Implementation of homomorphic filter made in Python, but something is wrong with it.
Source img: sourceimg
Using this may give me a filtered img like:filtered img Of course by changing the attributes I can make it lighter, but the excessive shading trend is obviously still exist.
Could anybody shed some light on this problem? Thx a lot!
Upvotes: 0
Views: 1490
Reputation: 114946
You don't say which version of SciPy that you are using, but I suspect you are hitting a bug that should have been fixed in SciPy version 1.0.1 (and 1.1.0). The bug causes the interpreter to crash when the input to the Gaussian filter function has data type numpy.float16
. If your input matrix fxy
has type numpy.uint8
(which is typical for an image), then numpy.log1p(fxy)
returns an array with data type numpy.float16
, e.g.:
In [29]: fxy = np.array([[0, 10], [128, 255]], dtype=np.uint8)
In [30]: logF = np.log1p(fxy)
In [31]: logF
Out[31]:
array([[0. , 2.398],
[4.86 , 5.547]], dtype=float16)
logF
has type numpy.float16
, and the crash occurs when you pass this to gaussian_filter1d
.
A quick fix is to convert fxy
to, say, 64 bit floating point before passing it to np.log1p
:
In [34]: logF = np.log1p(fxy.astype(np.float64))
In [35]: logF.dtype
Out[35]: dtype('float64')
Upvotes: 1