Add Phansalkar local threshold algorithm

[kleinerELM]

Applies Phansalkar local threshold to an array. Phansalkar is a modification of Sauvola technique to deal with low contrast images.

Description

I need the Phansalkar thresholding filter for low contrast SEM images. I found an implementation in ImageJ, which worked great. But I did not find any mentioning of that filter for scikit. I used threshold_sauvola() as a template.

Source paper:
Phansalskar N. et al. "Adaptive local thresholding for detection of nuclei in diversity stained cytology images.", International Conference on Communications and Signal Processing (ICCSP), pp. 218-220, 2011
DOI: 10.1109/ICCSP.2011.5739305

It should be noted, that the image has to be preprocessed using equalize_adapthist():

from skimage import data
from skimage.exposure import equalize_adapthist
image = data.moon()
image_eq = equalize_adapthist(image)
t_phansalkar = threshold_phansalkar(image_eq, window_size=15, k=0.25, p=2.0, q=10.0)
binary_image = image_eq > t_phansalkar

Checklist

• Docstrings for all functions
• Benchmark in ./benchmarks, if your changes aren't covered by an
  existing benchmark
• Gallery example in ./doc/examples (new features only)
• Unit tests
• Clean style in the spirit of PEP8

For reviewers

• Check that the PR title is short, concise, and will make sense 1 year
  later.
• Check that new functions are imported in corresponding __init__.py.
• Check that new features, API changes, and deprecations are mentioned in
  doc/release/release_dev.rst.

[pep8speaks]

Hello @kleinerELM! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2020-05-25 16:21:18 UTC

[lagru]

Thank you for the contribution @kleinerELM! This looks promising and from a first look at the paper and your description this might fit with the other thresholding functions. I'll give this a more thorough review later.

[rfezzani]

Thank you @kleinerELM for your contribution, I just left some comments ;)

[rfezzani]

Why not include this histogram equalization in the function since it seems that it is required for the method to work?

[kleinerELM]

I just followed the threshold_sauvola() fuction (and most of the other thresholding functions as well). They all return a threshold mask instead of a binary image. (Is there actually any reason for that?)
Since image_eq is necessary to create the correct binary image at the end, the function should best output a binary image directly, or what do you think?

I can see one advantage not to include the histogram equalization into the thresholding function:
If you want to find the best parameters window_size, k, r, p, q for a special image set, it is much faster if you have to calculate equalize_adapthist() only once, instead every time you change a parameter.
For example:

def threshold_phansalkar_parameter_variation( image ):
    #equalize only once
    image_eq = equalize_adapthist(image)
    
    # list of parameters
    k_list = [0.2,0.25,0.3,0.35]
    p_list = [1.5,2.0,2.5,3.0]

    fig, ax = plt.subplots(len(p_list),len(k_list), figsize=(15, 6) )
    fig.subplots_adjust(hspace = .5)
    ax = ax.ravel()

    for j in range(len(p_list)):
        for i in range(len(k_list)):
            # calculate binary mask of thresh_phansalkar() without doing equalize_adapthist() every time
            thresh_phansalkar = threshold_phansalkar(image_eq, window_size=15, k=k_list[i], p=p_list[j], q=10.0)
            binary_phansalkar = image_eq > thresh_phansalkar

            ax[i+j*len(k_list)].imshow(binary_phansalkar, cmap=plt.cm.gray)
            ax[i+j*len(k_list)].set_title( "k = " + str( k_list[i] )+ ", p = " + str( (p_list[j]) ) )
            ax[i+j*len(k_list)].set_axis_off()

    plt.show()

If you think that this is too rare a usecase, then your suggestion would definitely be more user-friendly.

[kleinerELM]

According to @jni the function should return a threshold map. Therefore, I think it is the best to leave the function as is. I documented the requirements of the function in the funtion itself and the example.

[rfezzani]

If the histogram equalization is integrated to the function, imax - imin can be replaced by image.ptp()

[kleinerELM]

Thank you for that suggestion.
I tried it for some images. And actually r will allways be 1, if the imput image is a result of the histogram equalization.
Therefore, I would suggest skipping the calculation and set r = 1.0 in that case.

What is the advantage of image.ptp()? I tried to benchmark the function and it seems image.ptp() was a little slower on my system than the current code. But that is negligible in comparison to the thresholding itself.

[lagru]

I think .ptp() it's just a shorthand. I find the current version (imax - imin) slightly more readable so I'd leave it as is personally but I don't mind either way.

[jni]

@kleinerELM btw this will need some tests, and I think it would be good to add to e.g. https://scikit-image.org/docs/dev/auto_examples/segmentation/plot_niblack_sauvola.html#sphx-glr-auto-examples-segmentation-plot-niblack-sauvola-py

It should definitely return a threshold map, not a binarised image.

[kleinerELM]

@kleinerELM btw this will need some tests, and I think it would be good to add to e.g. https://scikit-image.org/docs/dev/auto_examples/segmentation/plot_niblack_sauvola.html#sphx-glr-auto-examples-segmentation-plot-niblack-sauvola-py

It should definitely return a threshold map, not a binarised image.

Thank you for your patience. I added Phansalkar local threshold to the example and fixed the missing reference in skimage/filters/__init__.py .

While I was researching how to add the test functions I stumbled across the PR #2443 . It seems there was already a discussion including the Phansalkar threshold. Phansalkar was mistyped in that discussion, which is why I did not find that issue in the fist place.
Should I continue to work on the test suite or do we drop this PR? I personally would love to see this filter in Scikit.

[kleinerELM]

@jni I added test functions according to similar filter functions and they seem to work properly.
Please let me know if there's anything more I can do.

[kleinerELM]

Even if Phansalkar will not be accepted, I think it would be helpfull, if the most available filters are included in this function.

[kleinerELM]

I don't know if you would accept this as good style, but it works for thresholding algorithms which require different input variables.

[lagru]

Why not check for the actual function instead of a string comparison? E.g.

if func is threshold_local:
    ...

Also the parenthesis and extra space around the conditional expressions are a bit unconventional. I'm surprised that @pep8speaks isn't bugging you about it.

[kleinerELM]

Thanks for your advice, I changed it.
Actually, the parenthesis and extra space are a quirk of mine.

[lagru]

Am I correct in assuming that this line will always raise an AttributeError? At least that's what local testing tells me. Assuming you are the author, @sciunto do you still remember the intention behind this line?

I may be a good idea to use functools.wraps instead; wraps preserves the name of the function and makes the original one accessible through the __wrapped__ attribute.

[kleinerELM]

Can you make a code suggestion? I just followed the existing pattern here. I wanted to avoid breaking anything.

[lagru]

@kleinerELM I took the liberty to commit two small suggestions. You should be able to update your local branch with git pull origin threshold_phansalkar (in case you didn't know).