Non linear peak finding #84
Conversation
|
Is this necessary for this release? What about the updated PR for shm.py? I would much prefer to merge shm.py first. Let me know if the nonlinear has to be prioritized for some good reason. Great job, otherwise. |
| relative_peak_threshold, min_separation_angle) | ||
|
|
||
|
|
||
| def _nl_peak_finder(sphere_eval, seeds, fmin, gtol): | ||
| """Non-linear search for peaks from each seed |
There was a problem hiding this comment.
Let me see if I follow the logic of this: sphere_eval is a function that takes a Sphere as its input and returns a number. Here, we initialize the Sphere object with a single vector: theta, phi = x[0], x[1] and it tells us the value of this function in that location in spherical coordinates? In what follows, you get back a 'peak' of this function for every seed that we placed? That seems good. One things to think about is boundary conditions for theta and phi: what happens here when theta, phi go out of their domain (for example, above 2*pi for theta). The Sphere init doesn't seem to check for that, but maybe it should?
There was a problem hiding this comment.
You pretty much have the logic down. sphere_eval takes a sphere with N points and returns an array of length N. In this case we create a sphere with one point. The main peak finding function we have in dipy is called local_maxima. It is a discrete peak finding, we use it on line 140 to find the seeds for the non-linear search.
Boundary conditions are not an issue, I believe, because we eventually convert back to Cartesian coordinates and return unit vectors. ie theta=pi, phi={pi, 3pi} will both result in the same unit vector so I don't think we should should have any issues.
|
Do we have any peak-finding other than this, now that recspeed.peak_finding On Tue, Oct 9, 2012 at 11:41 AM, Eleftherios Garyfallidis <
|
|
The main peak finding function is dipy.reconst.recspeed.local_maxima it returns discrete peaks from a sphere, but sometimes you want the ability to find peaks with better precision than possible with a discrete sphere so in those cases this peak finder is useful. This peak finder makes many calls to |
| @@ -162,6 +162,50 @@ def remove_similar_vertices(cnp.ndarray[cnp.float_t, ndim=2, mode='strided'] ver | |||
| return unique_vertices[:count].copy() | |||
|
|
|||
|
|
|||
| @cython.boundscheck(False) | |||
| @cython.wraparound(False) | |||
| def search_descending(cnp.ndarray[cnp.float_t, ndim=1, mode='c'] a, | |||
There was a problem hiding this comment.
Is this faster than np.searchsorted?
There was a problem hiding this comment.
The short answer is yes. I tried using searchsorted but the overhead was too high, also it ends up looking ugly. The equivalent expression using search sorted is something like (len(a) - a[::-1].searchsorted(relative_threshold * a[0])) or 1. I figured having a dedicated function was cleaner and it allowed me to optimize it using cython.
There was a problem hiding this comment.
`searchsorted(.., side='left')`` should help
There was a problem hiding this comment.
I believe searchsorted uses side='left' as default, and I believe that's the side we need. The main issue is the len(a) - a[::-1] part adds a lot of overhead. I don't think numpy has a search descending only search ascending, ie searchsorted.
There was a problem hiding this comment.
Sounds like np.searchsorted(-x, -T).
There was a problem hiding this comment.
the biggest problem with this is that -x is an O(n) operation so you loose the benefit of doing a O(log(n)) search. It's all moot anyways for two reasons: 1, x tends to be small so this pretty fast no matter how you do it and 2, I like my little function and you're not going to take it away from me ;).
There was a problem hiding this comment.
OK, @MrBago knows my opinion on this one :)
|
This works beautifully. I coded up a version that optimized all the variables at the same time, but got exactly the same results without much (any?) speedup. When I initially tried this, my code broke due to the assumption that spheres store theta/phi as 1D arrays--I submitted a PR against your branch for that, @MrBago I'm +1 on merging. |
| # Remove small peaks | ||
| n = search_descending(values, relative_peak_threshold) | ||
| directions = vertices[:n] | ||
| # Remove peaks too close to each-other |
There was a problem hiding this comment.
On a purely aesthetic note, I prefer blank lines between logical sections, but I know this is contradictory to some other members' opinions, so I leave it to the author :)
|
Oh, I believe this warrants using one of the new GitHub emoticons: 👍 |
| self._config = {"sphere": sphere, | ||
| "relative_peak_threshold": relative_peak_threshold, | ||
| "min_separation_angle": min_separation_angle, | ||
| "fmin":fmin, "gtol":gtol} |
Ensure that sphere always has theta and phi coordinate packed as 1-D arrays.
|
@stefanv, when you have a sec try running this and see if you're still getting division by zero issues. |
|
@MrBago They're still there, unfortunately :/ |
…it is harder than I expected to support the fmin function as an argument to the NonLinearDirectionFinder. I've hard codded scipy.optimize.fmin as the optimization function for now and might re-visit this later.
We've been talking about this for a little while so, here it is!