Scikit-image: equalize_adapthist leaves artefacts when image size is not a multiple of tile size

[anntzer]

Python 3.5.2/Arch Linux/skimage 0.12.3

equalize_adapthist leaves "artefacts" when a dimension of the image is not a multiple of the kernel size:

from matplotlib import pyplot as plt
import numpy as np
from skimage.exposure import equalize_adapthist
img = (np.add.outer(np.sin(np.arange(126)), np.arange(128))
       + np.random.randint(3, size=(126, 128))).astype(int)
fig, axs = plt.subplots(1, 2)
axs[0].imshow(img, cmap="gray", interpolation="none")
axs[1].imshow(equalize_adapthist(img), cmap="gray", interpolation="none")
plt.show()

See the line at the top and the strip at the bottom of the second image.

Possibly related to #757?

[anntzer]

FWIW I now have an implementation of CLAHE which does not use any interpolation at all. For reasonably sized images such as this one, it is essentially instantaneous.

Would there be interest in merging it?

[soupault]

I'd definitely be glad to have a look.

[anntzer]

See https://github.com/anntzer/clahe. The implementation should be relatively transparent: we uniquify the values to make histogramming easier, and progressively update the histogram, and compute the CLAHE values. The only slightly tricky part is the contrast-limiting itself.
Edges are handled by padding with the border value.

[JDWarner]

@anntzer that looks rather elegant!

[anntzer]

I'll let you discuss the API changes :-)

Note that I use the "standard" definition of clip limit (>1, instead of <1 as is currently the case), see https://github.com/scikit-image/scikit-image/issues/2214#issuecomment-234841754.
For simplicity (mostly) and performance (less so), I only do one round of clipping, instead of repeating it until the contrast is actually limited by the specified value -- IMO there isn't much benefit in actually enforcing the exact clip limit.
Also I gave a try to "multiplicative" instead of "additive" redistribution of the clipped values (both are available), not an expert there but the original paper suggests that the latter (which is what's used in general) works better.

[b52]

What's the state on this one? This is a really annoying bug and it would be nice to see it fixed sooner rather than later, especially when there is a working code base available.

[soupault]

@anntzer I took a look into the details, and the implementation looks very promising. Could you manage to make it support multichannel images, and provide recommended default values for kernel_size and clip_limit?

[anntzer]

I'm not sure whether anything special is needed for multichannel images -- http://imagej.net/Enhance_Local_Contrast_(CLAHE) suggests that there is more to it than just processing each channel separately (which makes sense), but doesn't really go into the details of that; http://www.cs.unc.edu/Research/MIDAG/pubs/papers/Adaptive%20Histogram%20Equalization%20and%20Its%20Variations.pdf (which seems to be a good reference on the topic) doesn't go into multichannel images. (There seems to be more literature on the topic but I don't have the time to really go through it right now.)

It seems that a reasonable default for kernel_size is 1/8 of the image in both directions. I don't really see any standard default for clip_limit, and I'd think that defaulting to no contrast limitation is reasonable.

[soupault]

@anntzer alright!
Regarding multichannel images this appears to be one of the most popular approaches - https://stackoverflow.com/questions/24341114/simple-illumination-correction-in-images-opencv-c/24341809#24341809 , https://www.mathworks.com/help/images/ref/adapthisteq.html?requestedDomain=www.mathworks.com#bt56icn .
Would you like to work on these modifications?

[anntzer]

Upon further reflection, I believe that conversion to HSV/HSL/CIELAB/other color spaces really belongs to a separate preprocessing step, and that we should just document the possibilities (which could just be a three-line example -- conversion to color space, CLAHE, conversion back to original space -- in the docstring) rather than guess what's the best conversion for the user.

[soupault]

@anntzer I'd not agree. I believe that we do not provide any unified processing pipeline in, for example, CIE LAB/etc yet. Therefore, the algorithms have to be more or less self-sufficient.
What we can (and should, in my opinion) do here is to perform hist. equalization in the least desctructive (although, I'm not sure if AHE is possible in a non-destructive way) way (in RGB), and give to the user an option to select a better alternative (in CIE LAB).

You might look at the discussion on this topic in https://github.com/scikit-image/scikit-image/pull/2240 .

[jni]

@soupault although it would be very nice to have Lab conversion built into the function, I would certainly settle for a documentation fix in a PR if it'll move things along. imho the correct default thing to do would be to perform it in Lab — RGB CLAHE seems to me like it could cause terrible artifacts...

[anntzer]

I have no idea what "least destructive way" means. #2240 seems to suggest that the correct denoising space is dependent on the structure of the noise (which is not unexpected), so once again I'd rather not pick one for the user.

[soupault]

@anntzer In #2240 someone mentioned that, at the moment, the library functions are not meant to support as an input images in any possible colorspace. By default, we assume that the input is either grayscale or RGB[A] image.
In addition, different categories of image processing algorithms had been proven to perform better in _different_ colorspaces. From what I know: denoising - YCbCr/HSV/..., inpainting - LAB, HoG - RGB/LAB, etc. Of course, in a processing pipeline of an expert all colorspace conversions are usually performed manually, so we have to support the option to switch the conversion off. However, for beginners and advanced users, I think that an optimal processing out-of-the-box is better.

@jni

RGB CLAHE seems to me like it could cause terrible artifacts

Absolutely, that's why I'd like to avoid this case in as many situations as possible.

So, keeping in mind #2240, I'd propose:

def clahe(img, ..., convert2lab=True):
    ....
    img_conv = rgb2lab(img) if convert2lab else img
    img_clahe = _clahe(img_conv)
    img_out = lab2rgb(img_clahe) if convert2lab else img_clahe
    ....

We'll refactor this once metadata argument is added :D .

[anntzer]

Actually I realized I have another reason not to support multichannel images: the implementation I have currently also supports volumetric images (for which CLAHE is perfectly well defined -- and, empirically, works fine). So I'd rather interpret 3D arrays as volumes rather than RGB.

[jni]

@anntzer we love volumetric implementations in scikit-image! That's why we have the multichannel={True, False} keyword argument! =)

[anntzer]

I think I am willing to make a simple PR without multichannel support but if you insist on it I'll just leave it up to whoever is interested (feel free to reuse my code, of course).

[jni]

@anntzer I'm very happy to not have multichannel support! Please do submit it! I just meant that there is such a thing as supporting 3D, 2D+c, and 3D+c, and we have examples that do so. But 100% we would love a grayscale, nD CLAHE!

[anntzer]

Should it go in a separate function? I assume you'd want to keep the previous implementation, even if buggy, for backcompatibility...

[soupault]

@anntzer I think so.

[anntzer]

Hum, there actually seems to be some slight discrepancy between my Cython implementation and my pure-Python implementation... I'll need to investigate that first.

[anntzer]

It was just an off-by-1 error. Added some tests (that the two implementations, one on Python and one in Cython, do match).
@b52, would you be interested in making a PR out of my repo? I'm happy to walk you through my implementation, but I don't have time right now to write proper docs and sort out the machinery of how it'd interact with the currently existing version.

[b52]

I don't really have the time to create a PR out of something I only have a novice like understanding of.

[adamananos]

replying to the original post.

It appears that the problem is caused when the number of rows or number of columns computed from the kernel size provided to equalize_adapthist.clahe() results in an odd number. When the initial 'tile_size' - c_offset and r_offset are computed (first row or first column), this produces a non-integer (*.5). These values are rounded down in equalize_adapthist.interpolate() when they are converted to 'int'.

This can be corrected by adding a catch to see if the computed 'tile_size' value is divisible by 2 after the nr and nc are computed in equalize_adapthist.clahe() (see below)>>

nr = int(np.ceil(image.shape[0] / kernel_size[0]))
nc = int(np.ceil(image.shape[1] / kernel_size[1]))

add catches for the nr and nc to be divisible by 2
if nr%2 != 0:
nr -= 1
if nc%2 != 0:
nc -= 1

This appears to have been addressed before #757

[stefanv]

That looks like a regression; we should add a test to catch it. @adamananos Are you able to make a pull request to address the issue?

[RuRo]

Hi! Since you are going to reimplement/update adapthist, could you expose _adapthist.NR_OF_GREY as a parameter? If my understanding is correct, this determines the number of distinct gray colors used during equalization. I would like to use histogram normalization for medical images, which sometimes have 2**16 or even 2**32 distinct intensity values.

I am concerned, whether having NR_OF_GREY=2**14 hardcoded would result in the loss of details in the images.

[RuRo]

Unfortunately,

if nr%2 != 0:
    nr -= 1
if nc%2 != 0:
    nc -= 1

doesn't fix the issue, there are still some image sizes, where the black line appears.
If anyone knows, how to fix this or work around it, I would really appreciate it.

[YubinXie]

I am also seeing the weird lines in my image when use exposure.equalize_adapthist.
Now, I manually set the kernel size as submultiple of my image size to get rid of it. It would be nice to make it auto. This hidden artifacts make the method not reliable on some images and can be very dangerous.

[sciunto]

From the code written by the original poster in a jupyter notebook,

`from matplotlib import pyplot as plt
import numpy as np
%matplotlib inline
from skimage.exposure import equalize_adapthist
img = (np.add.outer(np.sin(np.arange(126)), np.arange(128))
+ np.random.randint(3, size=(126, 128))).astype(int)
fig, axs = plt.subplots(1, 2, figsize=(20, 20))
axs[0].imshow(img, cmap="gray", interpolation="none")
axs[1].imshow(equalize_adapthist(img), cmap="gray", interpolation="none")
plt.show()

from skimage import __version__
print(__version__)`

I obtain (sckimage 0.15.0).

Is the issue solved? Otherwise, can someone write a minimal code to reproduce?

[m-albert]

Hi!
I looked into the artefacts and I think there are two different types, which the example posted at the beginning of this issue nicely illustrates.
1) The light upper stripe. These artefacts can occur in the first tiles of every dimension and, as @adamananos says, come from having odd tile sizes. Dividing the odd tile size by two leads to non integer slices. The length of the slices are set here:
https://github.com/scikit-image/scikit-image/blob/2ffd9c823838a7975d1e16601cd724be530c31fa/skimage/exposure/_adapthist.py#L155-L156
and here:
https://github.com/scikit-image/scikit-image/blob/2ffd9c823838a7975d1e16601cd724be530c31fa/skimage/exposure/_adapthist.py#L169-L170
2) The lower dark stripe. This was addressed in #893 but then somehow reverted (as @stefanv
pointed out, can't seem to find where exactly). The tiled implementation clips the histogram only for tiles (of size kernel_size) naturally fitting into the shape of the image and leaves remaining pixels unprocessed. The tiles are defined here:
https://github.com/scikit-image/scikit-image/blob/2ffd9c823838a7975d1e16601cd724be530c31fa/skimage/exposure/_adapthist.py#L119-L123
In #893 this had been solved by padding the image to a multiple of the kernel_size before processing and then cropping the output image to the original shape.

If I'm not forgetting anything, both artefacts should be easy to eliminate:
1) Work on integer slices for the starting tiles np.ceil(kernel_size/2) and take care that the last tile doesn't overshoot the shape of the image.
2) Pad the image before processing and crop before returning the result, as done in #893

I added fixes for both artefacts in my generalisation of the existing CLAHE implementation to N dimensions (pull request here: #2761). Would there be interest in combining these fixes with the nd generalisation? I guess it's not super clean to mix things, but the code relevant for the fixes is different in the generalisation. Alternatively, I'm also happy to fix the artefacts for just the 2D version.

Here's the output of the code at the top in the fixed n-dimensional version:

Cheers!

[sciunto]

@jni, do you have a suggestion?

[m-albert]

Here's a quick side-by-side comparison between the skimage tiling approach and the window updating code of @anntzer at https://github.com/anntzer/clahe

from matplotlib import pylab as plt
import skimage
import clahe # https://github.com/anntzer/clahe
import time
from scipy import ndimage

np.random.seed(0)
imgb = np.random.random((10,10))
img = ndimage.zoom(imgb,100)
img = (img - img.min())/(img.max()-img.min())

ks = [10,20,50,100,200,500][:]
fig,axs = plt.subplots(len(ks),3,figsize = (6,12))
for ik,k in enumerate(ks):
    kernel_size = k
    clip_limit = 0.1

    axs[ik][0].imshow(img)
    axs[ik][0].set_title('k: %s' %k)
    tstart = time.time()
    res = skimage.exposure.equalize_adapthist((img*1000).astype(np.uint16),
                                              kernel_size=kernel_size,clip_limit=clip_limit)
    tend = time.time()
    axs[ik][1].imshow(res)
    axs[ik][1].set_title("{:3.2f}s".format(tend-tstart))
    tstart = time.time()
    res = clahe.clahe((img*1000).astype(np.uint16),
                      win_shape=[kernel_size]*2,clip_limit=1/clip_limit)
    tend = time.time()
    axs[ik][2].imshow(res)
    axs[ik][2].set_title("{:3.2f}s".format(tend-tstart))
plt.tight_layout()

comparison (visual and timing on a macbook)
skimage clahe vs github.com/anntzer/clahe
orig / skimage tiling / window updating

The window updating approach produces very smooth and pretty results. The tradeoff seems to be the processing time.

Here's another comparison in 3D using the branch from PR #2761:

from matplotlib import pylab as plt
import skimage
import clahe # https://github.com/anntzer/clahe
import time
from scipy import ndimage

np.random.seed(0)
imgb = np.random.random((5,5,5))
img = ndimage.zoom(imgb,40)
img = (img - img.min())/(img.max()-img.min())

ks = [10,20,50]
fig,axs = plt.subplots(len(ks),3,figsize = (6,2*len(ks)))
for ik,k in enumerate(ks):
    print(ik)
    kernel_size = k
    clip_limit = 0.1

    axs[ik][0].imshow(img[len(img)//2])
    axs[ik][0].set_title('k: %s' %k)
    tstart = time.time()
    res = skimage.exposure.equalize_adapthist((img*1000).astype(np.uint16),
                                              kernel_size=kernel_size,clip_limit=clip_limit)
    res = res[len(res)//2]
    tend = time.time()
    axs[ik][1].imshow(res)
    axs[ik][1].set_title("{:3.2f}s".format(tend-tstart))
    tstart = time.time()
    res = clahe.clahe((img*1000).astype(np.uint16),
                      win_shape=[kernel_size]*3,clip_limit=1/clip_limit)
    res = res[len(res)//2]
    tend = time.time()
    axs[ik][2].imshow(res)
    axs[ik][2].set_title("{:3.2f}s".format(tend-tstart))

fig.tight_layout(rect=[0, 0.05, 1.1, 0.95])
plt.tight_layout()

orig / skimage tiling / window updating

[rfezzani]

The described border artifacts in https://github.com/scikit-image/scikit-image/issues/2219#issuecomment-516526979 are fixed by #4349. May be opening an other issue or PR to implement @anntzer approach would be a good idea. Closing this one ;-) Thank you all.

[stefanv]

@rfezzani Please go ahead and open that issue for discussion, giving context.