Blob Detection ( DoG )

[vighneshbirodkar]

Example Code

from skimage.feature import get_blobs
from skimage import data
from matplotlib import pyplot as plt
import math

img = data.coins()
blobs = get_blobs(img)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
plt.imshow(img,cmap='gray')


for b in blobs:
    x,y = b[0],b[1]
    r = math.sqrt( b[2]/math.pi )
    c = plt.Circle((y,x),r,color='#ff0000',lw = 2,fill = False)
    ax.add_patch(c)

plt.show()

Results

This is my first PR to scikit-image. I would be glad if someone could go through it and suggest improvements

I haven't used Cython before.Could re writing _blob_overlap in Cython provide speed improvements ? The main bottle neck in this case is computing the Gaussian Filter.

@cdeil , @adonath Please share your thoughts.

LoG will be just minor additions to this. Once this is accepted I'll move onto that and the remaining ones.

[ahojnnes]

I have only had a quick look at the code, but it looks very nice and clean!

My recommendation API-wise: use the same convention as with skimage.feature.corner_*. There is certainly not only one blob detector. So we could call this one blob_dog, others would be called blob_xy.

[cdeil]

Thanks for implementing this!

Some comments:
The doctest for get_blobs is currently failing on travis-ci:
https://travis-ci.org/scikit-image/scikit-image/jobs/19771268#L2771

Clean up whitespace using the pep8 or autopep8 tools.
See complaints for blob.py starting here:
https://travis-ci.org/scikit-image/scikit-image/jobs/19771268#L2920

I don't have time to try this out today, but will soon.

[cdeil]

Missing comma, no?

[cdeil]

I've left some minor comments inline.

I think you should add some more tests, e.g. a test with two blobs in the image (currently there is only one test which has one blob).

[vighneshbirodkar]

@ahojnnes
The other blob detectors will be covered with the mode parameter. I could put in a function for each but a lot of things are common to most of the Blob Detectors mentioned here.
https://github.com/scikit-image/scikit-image/wiki/Requested-features

For eg : LoG can implemented with just two more lines after this
https://github.com/vighneshbirodkar/scikit-image/blob/master/skimage/feature/blob.py#L189

@cdeil
No I did not run the doctests locally. I will run them now and make changes. Thanks a lot for your help.

[adonath]

Of course I agree to use my code, it would be nice to have this functionality in scikit-image regularly. Thanks for working on this! The main credits should also go to: http://www.cs.utah.edu/~jfishbau/advimproc/project1/ I just adapted and modified the matlab code he provides on the bottom of the page. I will try to take a closer look at the code and comment within the next days.

[coveralls]

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

[vighneshbirodkar]

@adonath
I think it is better to have different functions for different algorithms, that allows different sets of default values

[cdeil]

Maybe split this formula into multiple statements?
(should improve speed, might increase memory)

Or at least improve formatting to make it more readable?

[cdeil]

I've left some very minor comments for typos.
Otherwise this looks pretty good to me!

Could one of the scikit-image devs please review it and comment on the API?

[cdeil]

What does the dog stand for?
Please mention this in the docstring and include a reference.

[vighneshbirodkar]

@cdeil
I have fixed all typos I could find.

@adonath
I followed this link to implement DoG
http://en.wikipedia.org/wiki/Blob_detection#The_difference_of_Gaussians_approach
Here they mention limit case implying dt to be very small

And while testing out I found that only keeping ds small I could get images similar to those obtained by gaussian_laplace

Regarding the point you made about Logarithmic scales. Consider this

>>> np.linspace(1,20,10)
array([  1.        ,   3.11111111,   5.22222222,   7.33333333,
         9.44444444,  11.55555556,  13.66666667,  15.77777778,
        17.88888889,  20.        ])
>>> np.logspace(math.log(1,10),math.log(20,10),10)
array([  1.        ,   1.39495079,   1.94588772,   2.71441762,
         3.78647901,   5.2819519 ,   7.368063  ,  10.27808533,
        14.33742329,  20.        ])

Using logarithmic scale results in fewer values in the higher ranges. If this were used, the size of larger blobs would be computed with less accuracy.
The user can increase max_sigma if he needs to find blobs with larger sizes

[cdeil]

Using logarithmic scale results in fewer values in the higher ranges. If this were used, the size of larger blobs would be computed with less accuracy.

Less absolute accuracy, yes.
But equal relative accuracy.

Imagine you have blobs in the size range 10 to 1000 and you want to to have ~ 10 % relative accuracy across that scale range. I think this is a common use case and is where you'd use logarithmic scales ... it's less common I think that you want to measure blobs of both size 10 and 1000 with a fixed absolute accuracy (of say 1).

[coveralls]

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

[stefanv]

Please note that we are now in the situation I have mentioned before, where
multiple images have been committed as part of this PR. In order to clean
it up, a rebase will be needed in which the commit with the wrong image is
removed in its entirety.

As for the image, using the original is better (we do not win anything by
converting it to PNG instead).

I am a tad concerned about the runtime here.

[vighneshbirodkar]

I'll send in a new PR with the JPG image, that's not a problem.

@stefanv
Regarding the runtime, even if no blobs are detected, multiple calls to gaussian_filter consume time.

[stefanv]

@vighneshbirodkar No, the problem is the growth in the repository size. We'll need to fix it in this PR.

[stefanv]

Ah, it looks like you haven't committed any images yet, so everything is fine.

[cdeil]

Since it's still not 100% clear which Hubble deep field example file will be included in the end, and there are ~ 30 commits in the history with typo / formatting changes, it might be simplest to, at the very end when this is ready, create a new branch off master and make one commit per file, no?
(This is what I would to here because I think it's simpler instead of removing and squashing commits from the history.)

[vighneshbirodkar]

@stefanv
What are you talking about ? There are 8 files in this PR, and one on them is an image

[stefanv]

@vighneshbirodkar My mistake. We'll have to rebase as I thought initially.

[vighneshbirodkar]

I'll wait to hear from the core developers.

Whose opinion should I wait for ?

[coveralls]

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

[cdeil]

Regarding the runtime, even if no blobs are detected, multiple calls to gaussian_filter consume time.

If I understand correctly, you currently run detection with:

blob_dog(gray_img, threshold=.5, min_sigma=0.1, max_sigma=10)

and the default parameters are

def blob_dog(image, min_sigma=1, max_sigma=20, num_sigma=50, delta=0.01,
             threshold=5.0, overlap=.5, log_scale=False):

I think you can use much fewer scales (and thus faster) if you use log_scale=True and get almost identical detections?

[vighneshbirodkar]

log_scale doesn't help, because the code will still compute 50*2 = 100 gaussians.

I put num_sigma=10 and got simlar results in ~1.4secs

[cdeil]

I would think log_scale=True would give better results here for a given num_sigma.
But OK ... as long as the detections visually roughly make sense (which they currently do) and the example runs in a second or two, it should be fine. Or @stefanv, would you prefer an even faster example?

[astrofrog]

@cdeil - I'd suggest asking on the astropy mailing list (broader than astropy-dev). Someone may already have one. Also, could you not merge this then substitute the image with AVM later on if needed?

[cdeil]

@stefanv Do you want to include the Hubble deep field file as it is now?
Or should we try and get one with AVM info (this defines the pixel to sky coordinate mapping) by asking on the astropy mailing list (which could mean a delay of a few days)?

Either way is fine with me ... AVM info is not really needed for scikit-image, which just uses pixel coordinates, it's just something that's nice to have for astronomers.

[stefanv]

@cdeil I think the PR is pretty much ready to go, but there is no sense of urgency to push the merge button. We can ask for the AVM info and hear what they say. But how would you include this meta-data? We cannot have FITS files in our data module, since pyfits is an optional dependency.

[cdeil]

AVM info can be included in JPEG or PNG files (see e.g. here).

[stefanv]

Perfect, let's get it then. I need to merge this one manually anyway to remove the old commits.

[vighneshbirodkar]

@stefanv
If you want I can create a new PR. If I need to include the AVM info, can you please send me the file @cdeil ?

[cdeil]

I've asked on the astropy mailing list:
http://mail.scipy.org/pipermail/astropy/2014-March/003085.html
Let's see if someone can simply point us to the best Hubble deep image file to include.

[cdeil]

I think we have a winner:
http://mail.scipy.org/pipermail/astropy/2014-March/003088.html

@stefanv That is if 570 kB is acceptable, it seems some of the other files in the repo are that large. If it's too large @vighneshbirodkar or I can browse around to find another one.

[vighneshbirodkar]

I'll crop out a region in the code

[adonath]

@vighneshbirodkar You already put a lot of work in this PR and it depends on how much effort you still want to put in it, but you might consider implementing the DoG blob detection using a Laplacian-Pyramid, which is much more efficient in terms of speed and memory, especially on larger images. The way it is implemented now it does not differ much (in terms of efficiency and speed) from the LoG approach, it could even be slower. You practically loose the main advantage of the DoG approach. This can be future work, but I think it should be done once the LoG blob detection is available. In my opinion the reason why you would choose the DoG over the LoG approach is because you would choose performance over accuracy of the detections.

[vighneshbirodkar]

@adonath
There are 5 requested algorithms on the wiki. I figured that the first 4 are similar, differing only in the way they compute each image in the scale space. So I thought once I had done any one of them, doing the others would be easier. I chose DoG because I had worked with it earlier. Once this gets accepted I am anyways going to implement the other 4. The last one ( MSER ) would require some Cython code I guess.

The SIFT algorithm as far as I know computes successive Gaussians with sigma being in the ratio of a constant k . I did try that, but I didn't get any results. The only times I got results were when took the difference between images using a small difference of sigma

[vighneshbirodkar]

@stefanv If I get a go ahead I'll use @cdeil 's image and create a new PR

[cdeil]

@vighneshbirodkar If you're making a new branch and PR anyways, it might be better to split the Hubble image new data set from the blob detect implementation.
This will make it easier to review and you'll probably get your code merged more quickly.
(That's actually my mistake, I should have opened a separate issue to suggest the new example image addition)

[vighneshbirodkar]

@cdeil , @stefanv
So I'll add the code and testcase is the next PR, and the example and Hubble image in the one after that ?

[vighneshbirodkar]

Resubmitted in #903