/usr/lib/python2.7/dist-packages/dipy/direction/probabilistic_direction_getter.py is in python-dipy 0.10.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 | """
Implementation of a probabilistic direction getter based on sampling from
discrete distribution (pmf) at each step of the tracking."""
import numpy as np
from .peaks import peak_directions, default_sphere
from dipy.reconst.shm import order_from_ncoef, sph_harm_lookup
from dipy.tracking.local.direction_getter import DirectionGetter
from dipy.tracking.local.interpolation import trilinear_interpolate4d
def _asarray(cython_memview):
# TODO: figure out the best way to get an array from a memory view.
# `np.array(view)` works, but is quite slow. Views are also "array_like",
# but using them as arrays seems to also be quite slow.
return np.fromiter(cython_memview, float)
class PmfGen(object):
pass
class SimplePmfGen(PmfGen):
def __init__(self, pmf_array):
if pmf_array.min() < 0:
raise ValueError("pmf should not have negative values")
self.pmf_array = pmf_array
def get_pmf(self, point):
return trilinear_interpolate4d(self.pmf_array, point)
class SHCoeffPmfGen(PmfGen):
def __init__(self, shcoeff, sphere, basis_type):
self.shcoeff = shcoeff
self.sphere = sphere
sh_order = order_from_ncoef(shcoeff.shape[-1])
try:
basis = sph_harm_lookup[basis_type]
except KeyError:
raise ValueError("%s is not a known basis type." % basis_type)
self._B, m, n = basis(sh_order, sphere.theta, sphere.phi)
def get_pmf(self, point):
coeff = trilinear_interpolate4d(self.shcoeff, point)
pmf = np.dot(self._B, coeff)
pmf.clip(0, out=pmf)
return pmf
class PeakDirectionGetter(DirectionGetter):
"""An abstract class for DirectionGetters that use the peak_directions
machinery."""
sphere = default_sphere
def __init__(self, sphere=None, **kwargs):
if sphere is not None:
self.sphere = sphere
self._pf_kwargs = kwargs
def _peak_directions(self, blob):
"""Gets directions using parameters provided at init.
Blob can be any function defined on ``self.sphere``, ie an ODF, PMF,
FOD.
"""
return peak_directions(blob, self.sphere, **self._pf_kwargs)[0]
class ProbabilisticDirectionGetter(PeakDirectionGetter):
"""Randomly samples direction of a sphere based on probability mass
function (pmf).
The main constructors for this class are current from_pmf and from_shcoeff.
The pmf gives the probability that each direction on the sphere should be
chosen as the next direction. To get the true pmf from the "raw pmf"
directions more than ``max_angle`` degrees from the incoming direction are
set to 0 and the result is normalized.
"""
@classmethod
def from_pmf(klass, pmf, max_angle, sphere, **kwargs):
"""Constructor for making a DirectionGetter from an array of Pmfs
Parameters
----------
pmf : array, 4d
The pmf to be used for tracking at each voxel.
max_angle : float, [0, 90]
The maximum allowed angle between incoming direction and new
direction.
sphere : Sphere
The set of directions to be used for tracking.
relative_peak_threshold : float in [0., 1.]
Used for extracting initial tracking directions. Passed to
peak_directions.
min_separation_angle : float in [0, 90]
Used for extracting initial tracking directions. Passed to
peak_directions.
See also
--------
dipy.direction.peaks.peak_directions
"""
pmf = np.asarray(pmf, dtype=float)
if pmf.ndim != 4:
raise ValueError("pmf should be a 4d array.")
if pmf.shape[3] != len(sphere.theta):
msg = ("The last dimension of pmf should match the number of "
"points in sphere.")
raise ValueError(msg)
pmf_gen = SimplePmfGen(pmf)
return klass(pmf_gen, max_angle, sphere, **kwargs)
@classmethod
def from_shcoeff(klass, shcoeff, max_angle, sphere, basis_type=None,
**kwargs):
"""Probabilistic direction getter from a distribution of directions
on the sphere.
Parameters
----------
shcoeff : array
The distribution of tracking directions at each voxel represented
as a function on the sphere using the real spherical harmonic
basis. For example the FOD of the Constrained Spherical
Deconvolution model can be used this way. This distribution will
be discretized using ``sphere`` and tracking directions will be
chosen from the vertices of ``sphere`` based on the distribution.
max_angle : float, [0, 90]
The maximum allowed angle between incoming direction and new
direction.
sphere : Sphere
The set of directions to be used for tracking.
basis_type : name of basis
The basis that ``shcoeff`` are associated with.
``dipy.reconst.shm.real_sym_sh_basis`` is used by default.
relative_peak_threshold : float in [0., 1.]
Used for extracting initial tracking directions. Passed to
peak_directions.
min_separation_angle : float in [0, 90]
Used for extracting initial tracking directions. Passed to
peak_directions.
See also
--------
dipy.direction.peaks.peak_directions
"""
pmf_gen = SHCoeffPmfGen(shcoeff, sphere, basis_type)
return klass(pmf_gen, max_angle, sphere, **kwargs)
def __init__(self, pmf_gen, max_angle, sphere=None, **kwargs):
"""Direction getter from a pmf generator.
Parameters
----------
pmf_gen : PmfGen
Used to get probability mass function for choosing tracking
directions.
max_angle : float, [0, 90]
The maximum allowed angle between incoming direction and new
direction.
sphere : Sphere
The set of directions to be used for tracking.
relative_peak_threshold : float in [0., 1.]
Used for extracting initial tracking directions. Passed to
peak_directions.
min_separation_angle : float in [0, 90]
Used for extracting initial tracking directions. Passed to
peak_directions.
See also
--------
dipy.direction.peaks.peak_directions
"""
PeakDirectionGetter.__init__(self, sphere, **kwargs)
self.pmf_gen = pmf_gen
# The vertices need to be in a contiguous array
self.vertices = self.sphere.vertices.copy()
cos_similarity = np.cos(np.deg2rad(max_angle))
self._set_adjacency_matrix(sphere, cos_similarity)
def _set_adjacency_matrix(self, sphere, cos_similarity):
"""Creates a dictionary where each key is a direction from sphere and
each value is a boolean array indicating which directions are less than
max_angle degrees from the key"""
matrix = np.dot(sphere.vertices, sphere.vertices.T)
matrix = abs(matrix) >= cos_similarity
keys = [tuple(v) for v in sphere.vertices]
adj_matrix = dict(zip(keys, matrix))
keys = [tuple(-v) for v in sphere.vertices]
adj_matrix.update(zip(keys, matrix))
self._adj_matrix = adj_matrix
def initial_direction(self, point):
"""Returns best directions at seed location to start tracking.
Parameters
----------
point : ndarray, shape (3,)
The point in an image at which to lookup tracking directions.
Returns
-------
directions : ndarray, shape (N, 3)
Possible tracking directions from point. ``N`` may be 0, all
directions should be unique.
"""
pmf = self.pmf_gen.get_pmf(point)
return self._peak_directions(pmf)
def get_direction(self, point, direction):
"""Samples a pmf to updates ``direction`` array with a new direction.
Parameters
----------
point : memory-view (or ndarray), shape (3,)
The point in an image at which to lookup tracking directions.
direction : memory-view (or ndarray), shape (3,)
Previous tracking direction.
Returns
-------
status : int
Returns 0 `direction` was updated with a new tracking direction, or
1 otherwise.
"""
# point and direction are passed in as cython memory views
pmf = self.pmf_gen.get_pmf(point)
cdf = (self._adj_matrix[tuple(direction)] * pmf).cumsum()
if cdf[-1] == 0:
return 1
random_sample = np.random.random() * cdf[-1]
idx = cdf.searchsorted(random_sample, 'right')
newdir = self.vertices[idx]
# Update direction and return 0 for error
if np.dot(newdir, _asarray(direction)) > 0:
direction[:] = newdir
else:
direction[:] = -newdir
return 0
class DeterministicMaximumDirectionGetter(ProbabilisticDirectionGetter):
"""Return direction of a sphere with the highest probability mass
function (pmf).
"""
def get_direction(self, point, direction):
"""Find direction with the highest pmf to updates ``direction`` array
with a new direction.
Parameters
----------
point : memory-view (or ndarray), shape (3,)
The point in an image at which to lookup tracking directions.
direction : memory-view (or ndarray), shape (3,)
Previous tracking direction.
Returns
-------
status : int
Returns 0 `direction` was updated with a new tracking direction, or
1 otherwise.
"""
# point and direction are passed in as cython memory views
pmf = self.pmf_gen.get_pmf(point)
cdf = self._adj_matrix[tuple(direction)] * pmf
idx = np.argmax(cdf)
if pmf[idx] == 0:
return 1
newdir = self.vertices[idx]
# Update direction and return 0 for error
if np.dot(newdir, _asarray(direction)) > 0:
direction[:] = newdir
else:
direction[:] = -newdir
return 0
|