Blob Detection ( DoG ) - Resubmitting

[vighneshbirodkar]

Resubmitting #900

As @adonath has mentioned, the last pull request didn't provide significant speed improvements over the Laplacian Of Gaussians (LoG) approach. This computes half the number of Gaussians the last method and is significantly faster than the LoG approach. But the trade-off is accuracy. The esitmated size of the blobs is less accurate than the previous method. LoG approach provides more accurate size estimation ( which I will implement later )

I haven't added an example here. I'll add it in my next PR after this is merged, which will include a example with the Hubble Deep Field image.

Examples

[coveralls]

Coverage remained the same when pulling 59fa47c on vighneshbirodkar:dog into d652a9c on scikit-image:master.

[coveralls]

Coverage remained the same when pulling 59fa47c on vighneshbirodkar:dog into d652a9c on scikit-image:master.

[coveralls]

Coverage remained the same when pulling 5a988fd on vighneshbirodkar:dog into d652a9c on scikit-image:master.

[jni]

I think this function can and should be made nd and take a connectivity parameter (so you can use 6 neighbors, 18, or 26; I prefer to refer to these as connectivity = 1, 2, 3), by using scipy.ndimage.generate_binary_structure. Then it might be useful as a public function.

[vighneshbirodkar]

How should we keep it in the API ? skimage.feature.get_local_maxima_3d ?

[vighneshbirodkar]

@jni
I made the changes . How should I expose it ? skimage.feature.get_local_maxima_3d ?

[coveralls]

Coverage remained the same when pulling e449ce9 on vighneshbirodkar:dog into d652a9c on scikit-image:master.

[jni]

@vighneshbirodkar, thanks, additional changes:

• change array to more commonly used im, img or ar. (array can be confused with the class np.array.)
• use ar.ndim instead of 3 as the rank, and rename function to just get_local_maxima, which now works in nd.
• change docstring under connectivity to specify that number of neighbors is an example for 3D arrays.
• I think exposing this under feature is not bad, but not great either. I might favor morphology. @stefanv what do you think? Pigeonholes are hard... =)
• Maybe add a doctest example showing the output for a 2D, 5x5 image with two local maxima.

[ahojnnes]

Why not incorporate this into skimage.feature.peak_local_max?

[vighneshbirodkar]

@ahojnnes I think peak_local_max is doing a lot of extra things which are unnecessary for the specific operation of blob detection. And it's also doing a lot of operations assuming that the input is a 2D array

[vighneshbirodkar]

Should I leave it here and import it in skimage/morphology/__init__.py or move it elsewhere ?

[tonysyu]

Should I leave it here and import it in skimage/morphology/init.py or move it elsewhere ?

The function should only be made available from one place, otherwise it gets confusing. I would favor keeping it in feature rather than moving it to morphology: I think of morphology functions as manipulation and feature functions as detection (there are probably contradictions on both sides, but whatever).

BTW, thanks for all your work on this PR! I've been watching from the sidelines. I know it can be a lot of work the first time around, but the review process really improves the quality of the code base and makes it much easier for everyone who reads it weeks, months, or years down the road.

[vighneshbirodkar]

So for now I have kept it as skimage.feature.get_local_maxima. If we were to move it to morphology the best place for it to go would be inside skimage/morphology/peak.py

@tonysyu , @jni , @stefanv your thoughts ?

[stefanv]

@ahojnnes Do you have any opinion about the peak finder?

[vighneshbirodkar]

I was also thinking of implementing other blob detection algorithms mentioned in the wiki. Should I add them to this PR ? Or wait for this one to get merged and add them in the next ?

[ahojnnes]

Sorry for asking again, but what does this function provide that skimage.feature.peak_local_max does not provide?

[vighneshbirodkar]

skimage.feature.peak_local_max assumes that input is an image, where as get_local_maxima can work with any array ( a 3d array in case of blob detections )

[ahojnnes]

In my opinion we should add the 3D functionality to peak_local_max... There is no reason to separate this functionality in two functions?

[vighneshbirodkar]

peak_local_max does a lot of other things than just find peaks, and we don't want it just for the 3D case, as @jni said, if the function operates on an ND array, it's useful as a public function.

[stefanv]

I guess the question is whether it is worth duplicating functionality here.
Would it be possible to update the existing peak finder for N-d, or would
that be prohibitively expensive?

[vighneshbirodkar]

I am guessing that is the case, especially handing the case where peak_local_max excludes the border of the input array

[ahojnnes]

Sorry, but I still don't understand why peak_local_max would not work with N-D arrays... I just tested:

In [1]: from skimage.feature import peak_local_max

In [2]: import numpy as np

In [3]: a = np.zeros((20, 20, 20))

In [4]: a[10, 10, 10] = 1

In [5]: peak_local_max(a, min_distance=1)
Out[5]: array([[10, 10, 10]])

Though, when looking at the code, there might be a bug in the logic of exclude_border, shouldn't we set the values to -INF in lines 136, 137?

[vighneshbirodkar]

@ahojnnes You are right, I switched to peak_local_max, with the right parameters, results are identical

[coveralls]

Coverage remained the same when pulling 7986714 on vighneshbirodkar:dog into d652a9c on scikit-image:master.

[coveralls]

Coverage remained the same when pulling 7986714 on vighneshbirodkar:dog into d652a9c on scikit-image:master.

[jni]

@everyone, if peak_local_max works in nd, the documentation of that
function should be updated. (Output doc at least assumes 2D.)

On Sunday, March 9, 2014, Coveralls notifications@github.com wrote:

[image: Coverage Status] https://coveralls.io/builds/582532

Coverage remained the same when pulling 7986714
7986714
on vighneshbirodkar:dog into d652a9c
d652a9c
on scikit-image:master.

Reply to this email directly or view it on GitHubhttps://github.com//pull/903#issuecomment-37131525
.

[stefanv]

@jni Please file an issue to keep track.

[stefanv]

I think we are ready to merge. Please replace these two lines with entries for yourself and @adonath in CONTRIBUTORS.txt.

[vighneshbirodkar]

@stefanv I will implement the remaining blob detection algorithms in subsequent PRs and then add an example

[coveralls]

Coverage remained the same when pulling 0f44e6f on vighneshbirodkar:dog into d652a9c on scikit-image:master.

[stefanv]

Thanks, @vighneshbirodkar! Fantastic work.

[cdeil]

Thank you!

Please cc me if you make a new PR with the Hubble deep field example later.

[ahojnnes]

@jni See #906.

[ankit-maverick]

@stefanv @vighneshbirodkar : I went through the conversation in the previous PR about returning sigma and I think it should be returned as well as the fourth column. Its true that the area is more generic while sigma is specific to specific blob detectors, but then again I guess we are having different blob_* functions for different blob detectors and it's better to preserve information instead of throwing it :).

[vighneshbirodkar]

@ankit-maverick
The area is just a function of sigma. Also do you think any case where sigma would acutally be useful to any user who is doing blob detection ?

[ankit-maverick]

@vighneshbirodkar : Yes, the area and sigma have one to one correspondence and I agree it is pretty direct. A case where a user might want sigma is for interpretation in sigma-domain(like a function in space or time can be equivalently representation in frequency domain but the kind of interpretation required may depend on the task). For example, the detector in SIFT(https://en.wikipedia.org/wiki/Scale-invariant_feature_transform) uses DoG!! and if you go through the SIFT paper(http://www.cs.ubc.ca/~lowe/papers/ijcv04.pdf), the scale-space(that introduces the scale invariance) discussion proceeds in terms of sigma and octaves because it is easier to do so instead of discussing and formulating using area.

[vighneshbirodkar]

@ankit-maverick
I don't think this function should be used as a part of something like SIFT. Although the DoG approach is inspired from SIFT it does not do the exact same thing as SIFT. If you notice, the code here multiplies the difference with sigma to provide scale invariance while SIFT does not do so. Therefore the using this as a part inside SIFT will produce different results

[ankit-maverick]

@vighneshbirodkar

I don't think this function should be used as a part of something like SIFT.

True.

Although the DoG approach is inspired from SIFT,

If I am not wrong, it's the other way round :)

it does not do the exact same thing as SIFT.....

I agree that this DoG blob detector is not intended to be used inside SIFT, but my point was that the sigma values provide a better interpretation(intuition?) as to what is happening in the scale space. Area, for sure is a better metric when this blob detector is used for application purposes.

[vighneshbirodkar]

@ankit-maverick
It can be done, but returning an array where one value is a function of the other seems a bit odd to me, we'll hear what the others have to say

[cdeil]

I like how regionprops returns the info ... with a simple array its easy to confuse which column is which.
In this case area or sigma could be stored and the other could be a property.

But before someone from the scikit-image team commented that they generally prefer simple arrays, and in this case there's much fewer properties of the detected objects, so an array is not that bad.
If we stay with the array as return value, I have no preference for area and / or sigma.

[stefanv]

As I might have mentioned during the last discussion, my preference is to
return sigma. I find the area less intuitive, ironically. I presume each
blob detector should return the "model parameters" involved, which would
typically not be area.

[vighneshbirodkar]

Working on it, I'll submit a PR with the change, I've already implemented LoG algorithm, so I'll include that as well.I'll remove area and return just (x,y,sigma)