Fix zero division error in SIFT #7069
Conversation
skimage/feature/sift.py
Outdated
| @@ -429,7 +429,7 @@ def _find_localize_evaluate(self, dogspace, img_shape): | |||
|
|
|||
| def _fit(self, h): | |||
| """Refine the position of the peak by fitting it to a parabola""" | |||
| return (h[0] - h[2]) / (2 * (h[0] + h[2] - 2 * h[1])) | |||
| return (h[0] - h[2]) / (2 * ((h[0] - h[1]) + (h[2] - h[1]))) | |||
There was a problem hiding this comment.
Arithmetically, it should be exactly the same, but this change prevents "zero" most of the time if not all because h[0] > = h[1] and h[2] >= h[1] are guaranteed and either of them should not hold equality.
There was a problem hiding this comment.
Makes sense. Thanks for reporting. Would you mind adding a comment above this line?
There was a problem hiding this comment.
Sure, I can add a comment.
There was a problem hiding this comment.
Thanks, I'll wait for your new commit.
skimage/feature/sift.py
Outdated
| @@ -429,7 +429,7 @@ def _find_localize_evaluate(self, dogspace, img_shape): | |||
|
|
|||
| def _fit(self, h): | |||
| """Refine the position of the peak by fitting it to a parabola""" | |||
| return (h[0] - h[2]) / (2 * (h[0] + h[2] - 2 * h[1])) | |||
| return (h[0] - h[2]) / (2 * ((h[0] - h[1]) + (h[2] - h[1]))) | |||
There was a problem hiding this comment.
Makes sense. Thanks for reporting. Would you mind adding a comment above this line?
skimage/feature/sift.py
Outdated
| @@ -519,6 +519,8 @@ def _compute_orientation(self, positions_oct, scales_oct, sigmas_oct, | |||
| # result | |||
| ori = ((m + self._fit(hist[neigh]) + 0.5) | |||
| * 2 * np.pi / self.n_bins) | |||
| if not np.isfinite(ori): | |||
| continue | |||
There was a problem hiding this comment.
Shouldn't the exception be handled explicitly (i.e., with a try statement and except ZeroDivisionError)? And within the _fit function?
There was a problem hiding this comment.
Let's say if we handle it by exception within _fit. What do you want to return by _fit when we catch the exception?
There was a problem hiding this comment.
I admit I would let it error, with a helpful message though: Maybe suggesting that the user pass h with a different precision, to circumvent the ZeroDivisionError?
There was a problem hiding this comment.
Hmm, sorry, I guess h is not really 'passed by the user' at this point; it's usually an intermediate result in another private function... Maybe try returning h[0] - h[2]) / (2 * ((h[0] - h[1]) + (h[2] - h[1])) as you wrote, and if this fails return h[0] - h[2]) / (2 * (h[0] + h[2] - 2 * h[1]) as originally, so we have better chances overall...? And handle the exception in _compute_orientation instead?
There was a problem hiding this comment.
Sorry for the late reply. I was too busy to take a look on it.
Maybe try returning h[0] - h[2]) / (2 * ((h[0] - h[1]) + (h[2] - h[1])) as you wrote, and if this fails return h[0] - h[2]) / (2 * (h[0] + h[2] - 2 * h[1])
I don't think it increases the chance. But, in any case, I guess it could be h[0] == h[1] == h[2] because
In [2]: ndi.maximum_filter([1., 1., 1.], [3], mode='wrap')
Out[2]: array([1., 1., 1.])So, I think ori could be nan and inf, and isfinite can handle both cases.
And handle the exception in _compute_orientation instead?
Actually, my PR description was inaccurate. The zero division was actually happening, but the error happened at this line with ValueError: math domain error because angle became inf. (Maybe better to replace cos and sin to numpy's function.)
So, unless we cast each element of h to python built-in type, we can't catch it as ZeroDivisionError because h is a numpy array.
By the way, I have two more PRs open for fixing bugs of SIFT (#7065 and #6862). So, it would be great if you can review them too.
There was a problem hiding this comment.
Actually, my PR description was inaccurate.
No worries; feel free to edit your PR's title.
Actually, my PR description was inaccurate. The zero division was actually happening, but the error happened at this line with ValueError: math domain error because angle became inf.
Thanks for investigating, @hayatoikoma. So, I think that your addition line 522 would be a good enough fix... How about reverting your change line 432?
(Maybe better to replace cos and sin to numpy's function.)
I agree.
There was a problem hiding this comment.
Yes, it should handle both, but why do you want to revert line 432?
There was a problem hiding this comment.
Can you please add a test? Maybe you could write a toy version of your use case? So it breaks with current main, and passes with your proposed changes.
Yes, it should handle both
Sorry, I don't get it: both what?
but why do you want to revert line 432?
Let's not reverse the burden of proof... 😉 Unless you are fixing two different things and I didn't get it. I was under the impression that the change line 522 was enough because it handles the case when the angle becomes inf or nan. If you are changing line 432 as well, please add a commend and/or test for it. Thanks for following through!
There was a problem hiding this comment.
Thinking about it again, I thought I should just handle when the denominator becomes zero. Then, it should handle the cases I reported. What do you think?
Description
I encountered a math domain error at
cosin_rotatebecause_fitreturnedinfwith float32 due to numerical precision.In the failed case,
h[0]was slightly smaller thanh[1]buth[1]andh[2]is exactly the same. Arithmetically, the divisor should not become zero, but it became exactly zero due to numerical precision.Release note