scikit-image

1
#cython: cdivision=True
2
#cython: nonecheck=False
3
"""
4
Cython implementation of Dijkstra's minimum cost path algorithm,
5
for use with data on a n-dimensional lattice.
6
"""
7
import cython
8
import numpy as np
9
from . import heap
10
from .._shared.utils import warn
11

12
cimport numpy as cnp
13
from . cimport heap
14

15
cnp.import_array()
16

17
OFFSET_D = np.int8
18
OFFSETS_INDEX_D = np.int16
19
EDGE_D = np.int8
20
INDEX_D = np.intp
21
FLOAT_D = np.float64
22

23

24
@cython.boundscheck(False)
25
@cython.wraparound(False)
26
def _get_edge_map(shape):
27
    """Return an array with edge points/lines/planes/hyperplanes marked.
28

29
    Given a shape (of length n), return an edge_map array with a shape of
30
    original_shape + (n,), where, for each dimension, edge_map[...,dim] will
31
    have zeros at indices not along an edge in that dimension, -1s at indices
32
    along the lower boundary, and +1s on the upper boundary.
33

34
    This allows one to, given an nd index, calculate not only if the index is
35
    at the edge of the array, but if so, which edge(s) it lies along.
36

37
    """
38
    d = len(shape)
39
    edges = np.zeros(shape+(d,), order='F', dtype=EDGE_D)
40
    for i in range(d):
41
        slices = [slice(None)] * (d+1)
42
        slices[d] = i
43
        slices[i] = 0
44
        edges[tuple(slices)] = -1
45
        slices[i] = -1
46
        edges[tuple(slices)] = 1
47
    return edges
48

49

50
@cython.boundscheck(False)
51
@cython.wraparound(False)
52
def _offset_edge_map(shape, offsets):
53
    """Return an array with positions marked where offsets will step
54
    out of bounds.
55

56
    Given a shape (of length n) and a list of n-d offsets, return a two arrays
57
    of (n,) + shape: pos_edge_map and neg_edge_map.
58
    For each dimension xxx_edge_map[dim, ...] has zeros at indices at which
59
    none of the given offsets (in that dimension) of the given sign (positive
60
    or negative, respectively) will step out of bounds. If the value is
61
    nonzero, it gives the largest offset (in terms of absolute value) that
62
    will step out of bounds in that direction.
63

64
    An example will be explanatory:
65
    >>> offsets = [[-2,0], [1,1], [0,2]]
66
    >>> pos_edge_map, neg_edge_map = _offset_edge_map((4,4), offsets)
67
    >>> neg_edge_map[0]
68
    array([[-1, -1, -1, -1],
69
          [-2, -2, -2, -2],
70
          [ 0,  0,  0,  0],
71
          [ 0,  0,  0,  0]], dtype=int8)
72

73
    >>> pos_edge_map[1]
74
    array([[0, 0, 2, 1],
75
          [0, 0, 2, 1],
76
          [0, 0, 2, 1],
77
          [0, 0, 2, 1]], dtype=int8)
78

79
    """
80
    indices = np.indices(shape)  # indices.shape = (n,)+shape
81

82
    #get the distance from each index to the upper or lower edge in each dim
83
    pos_edges = (shape - indices.T).T
84
    neg_edges = -1 - indices
85
    # now set the distances to zero if none of the given offsets could reach
86
    offsets = np.asarray(offsets)
87
    maxes = offsets.max(axis=0)
88
    mins = offsets.min(axis=0)
89
    for pos, neg, mx, mn in zip(pos_edges, neg_edges, maxes, mins):
90
        pos[pos > mx] = 0
91
        neg[neg < mn] = 0
92
    return pos_edges.astype(EDGE_D), neg_edges.astype(EDGE_D)
93

94

95
@cython.boundscheck(False)
96
@cython.wraparound(False)
97
def make_offsets(d, fully_connected):
98
    """Make a list of offsets from a center point defining a n-dim
99
    neighborhood.
100

101
    Parameters
102
    ----------
103
    d : int
104
        dimension of the offsets to produce
105
    fully_connected : bool
106
        whether the neighborhood should be singly- of fully-connected
107

108
    Returns
109
    -------
110
    offsets : list of tuples of length `d`
111

112
    Examples
113
    --------
114

115
    The singly-connected 2-d neighborhood is four offsets:
116

117
    >>> make_offsets(2, False)
118
    [(-1,0), (1,0), (0,-1), (0,1)]
119

120
    While the fully-connected 2-d neighborhood is the full cartesian product
121
    of {-1, 0, 1} (less the origin (0,0)).
122

123
    """
124
    if fully_connected:
125
        mask = np.ones([3]*d, dtype=np.uint8)
126
        mask[tuple([1]*d)] = 0
127
    else:
128
        mask = np.zeros([3]*d, dtype=np.uint8)
129
        for i in range(d):
130
            indices = [1]*d
131
            indices[i] = (0, -1)
132
            mask[tuple(indices)] = 1
133
    offsets = []
134
    for indices, value in np.ndenumerate(mask):
135
        if value == 1:
136
            indices = np.array(indices) - 1
137
            offsets.append(indices)
138
    return offsets
139

140

141
@cython.boundscheck(True)
142
@cython.wraparound(True)
143
def _unravel_index_fortran(flat_indices, shape):
144
    """_unravel_index_fortran(flat_indices, shape)
145

146
    Given a flat index into an n-d fortran-strided array, return an
147
    index tuple.
148

149
    """
150
    strides = np.multiply.accumulate([1] + list(shape[:-1]))
151
    indices = [tuple((idx // strides) % shape) for idx in flat_indices]
152
    return indices
153

154

155
@cython.boundscheck(True)
156
@cython.wraparound(True)
157
def _ravel_index_fortran(indices, shape):
158
    """_ravel_index_fortran(flat_indices, shape)
159

160
    Given an index tuple into an n-d fortran-strided array, return a
161
    flat index.
162

163
    """
164
    strides = np.multiply.accumulate([1] + list(shape[:-1]))
165
    flat_indices = [np.sum(strides * idx) for idx in indices]
166
    return flat_indices
167

168

169
@cython.boundscheck(False)
170
@cython.wraparound(False)
171
def _normalize_indices(indices, shape):
172
    """_normalize_indices(indices, shape)
173

174
    Make all indices positive. If an index is out-of-bounds, return None.
175

176
    """
177
    new_indices = []
178
    for index in indices:
179
        if len(index) != len(shape):
180
            return None
181
        new_index = []
182
        for i, s in zip(index, shape):
183
            i = int(i)
184
            if i < 0:
185
                i = s + i
186
            if not (0 <= i < s):
187
                return None
188
            new_index.append(i)
189
        new_indices.append(new_index)
190
    return new_indices
191

192

193
@cython.boundscheck(True)
194
@cython.wraparound(True)
195
def _reverse(arr):
196
    """Reverse index an array safely, with bounds/wraparound checks on.
197
    """
198
    return arr[::-1]
199

200

201
@cython.boundscheck(False)
202
@cython.wraparound(False)
203
cdef class MCP:
204
    """MCP(costs, offsets=None, fully_connected=True, sampling=None)
205

206
    A class for finding the minimum cost path through a given n-d costs array.
207

208
    Given an n-d costs array, this class can be used to find the minimum-cost
209
    path through that array from any set of points to any other set of points.
210
    Basic usage is to initialize the class and call find_costs() with a one
211
    or more starting indices (and an optional list of end indices). After
212
    that, call traceback() one or more times to find the path from any given
213
    end-position to the closest starting index. New paths through the same
214
    costs array can be found by calling find_costs() repeatedly.
215

216
    The cost of a path is calculated simply as the sum of the values of the
217
    `costs` array at each point on the path. The class MCP_Geometric, on the
218
    other hand, accounts for the fact that diagonal vs. axial moves are of
219
    different lengths, and weights the path cost accordingly.
220

221
    Array elements with infinite or negative costs will simply be ignored, as
222
    will paths whose cumulative cost overflows to infinite.
223

224
    Parameters
225
    ----------
226
    costs : ndarray
227
    offsets : iterable, optional
228
        A list of offset tuples: each offset specifies a valid move from a
229
        given n-d position.
230
        If not provided, offsets corresponding to a singly- or fully-connected
231
        n-d neighborhood will be constructed with make_offsets(), using the
232
        `fully_connected` parameter value.
233
    fully_connected : bool, optional
234
        If no `offsets` are provided, this determines the connectivity of the
235
        generated neighborhood. If true, the path may go along diagonals
236
        between elements of the `costs` array; otherwise only axial moves are
237
        permitted.
238
    sampling : tuple, optional
239
        For each dimension, specifies the distance between two cells/voxels.
240
        If not given or None, the distance is assumed unit.
241

242
    Attributes
243
    ----------
244
    offsets : ndarray
245
        Equivalent to the `offsets` provided to the constructor, or if none
246
        were so provided, the offsets created for the requested n-d
247
        neighborhood. These are useful for interpreting the `traceback` array
248
        returned by the find_costs() method.
249

250
    """
251

252
    def __init__(self, costs, offsets=None, fully_connected=True,
253
                 sampling=None):
254
        """__init__(costs, offsets=None, fully_connected=True, sampling=None)
255

256
        See class documentation.
257
        """
258
        costs = np.asarray(costs)
259
        if not np.can_cast(costs.dtype, FLOAT_D):
260
            raise TypeError('cannot cast costs array to ' + str(FLOAT_D))
261

262
        # Check sampling
263
        if sampling is None:
264
            sampling = np.array([1.0 for s in costs.shape], FLOAT_D)
265
        elif isinstance(sampling, (list, tuple)):
266
            sampling = np.array(sampling, FLOAT_D)
267
            if sampling.ndim != 1 or len(sampling) != costs.ndim:
268
                raise ValueError('Need one sampling element per dimension.')
269
        else:
270
            raise ValueError('Invalid type for sampling: %r.' % type(sampling))
271

272
        # We use flat, fortran-style indexing here (could use C-style,
273
        # but this is my code and I like fortran-style! Also, it's
274
        # faster when working with image arrays, which are often
275
        # already fortran-strided.)
276
        try:
277
            self.flat_costs = costs.astype(FLOAT_D, copy=False).ravel('F')
278
        except TypeError:
279
            self.flat_costs = costs.astype(FLOAT_D).flatten('F')
280
            warn('Upgrading NumPy should decrease memory usage and increase'
281
                 ' speed.')
282
        size = self.flat_costs.shape[0]
283
        self.flat_cumulative_costs = np.empty(size, dtype=FLOAT_D)
284
        self.dim = len(costs.shape)
285
        self.costs_shape = costs.shape
286
        self.costs_heap = heap.FastUpdateBinaryHeap(initial_capacity=128,
287
                                                    max_reference=size-1)
288

289
        # This array stores, for each point, the index into the offset
290
        # array (see below) that leads to that point from the
291
        # predecessor point.
292
        self.traceback_offsets = np.empty(size, dtype=OFFSETS_INDEX_D)
293

294
        # The offsets are a list of relative offsets from a central
295
        # point to each point in the relevant neighborhood. (e.g. (-1,
296
        # 0) might be a 2d offset).
297
        # These offsets are raveled to provide flat, 1d offsets that can be
298
        # used in the same way for flat indices to move to neighboring points.
299
        if offsets is None:
300
            offsets = make_offsets(self.dim, fully_connected)
301
        self.offsets = np.array(offsets, dtype=OFFSET_D)
302
        self.flat_offsets = np.array(
303
            _ravel_index_fortran(self.offsets, self.costs_shape),
304
            dtype=INDEX_D)
305

306
        # Instead of unraveling each index during the pathfinding algorithm, we
307
        # will use a pre-computed "edge map" that specifies for each dimension
308
        # whether a given index is on a lower or upper boundary (or none at
309
        # all). Flatten this map to get something that can be indexed as by the
310
        # same flat indices as elsewhere.
311
        # The edge map stores more than a boolean "on some edge" flag so as to
312
        # allow us to examine the non-out-of-bounds neighbors for a given edge
313
        # point while excluding the neighbors which are outside the array.
314
        pos, neg = _offset_edge_map(costs.shape, self.offsets)
315
        self.flat_pos_edge_map = pos.reshape((self.dim, size), order='F')
316
        self.flat_neg_edge_map = neg.reshape((self.dim, size), order='F')
317

318

319
        # The offset lengths are the distances traveled along each offset
320
        self.offset_lengths = np.sqrt(np.sum((sampling * self.offsets)**2,
321
                                      axis=1)).astype(FLOAT_D)
322
        self.dirty = 0
323
        self.use_start_cost = 1
324

325

326
    def _reset(self):
327
        """_reset()
328
        Clears paths found by find_costs().
329
        """
330

331
        cdef INDEX_T start
332

333
        self.costs_heap.reset()
334
        self.traceback_offsets[...] = -2  # -2 is not reached, -1 is start
335
        self.flat_cumulative_costs[...] = np.inf
336
        self.dirty = 0
337

338
        # Get starts and ends
339
        # We do not pass them in as arguments for backwards compat
340
        starts, _ = self._starts, self._ends
341

342
        # push each start point into the heap. Note that we use flat indexing!
343
        for start in _ravel_index_fortran(starts, self.costs_shape):
344
            self.traceback_offsets[start] = -1
345
            if self.use_start_cost:
346
                self.costs_heap.push_fast(self.flat_costs[start], start)
347
            else:
348
                self.costs_heap.push_fast(0, start)
349

350

351
    cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost,
352
                              FLOAT_T new_cost, FLOAT_T offset_length):
353
        """ float _travel_cost(float old_cost, float new_cost,
354
                               float offset_length)
355
        The travel cost for going from the current node to the next.
356
        Default is simply the cost of the next node.
357
        """
358
        return new_cost
359

360

361
    cpdef int goal_reached(self, INDEX_T index, FLOAT_T cumcost):
362
        """ int goal_reached(int index, float cumcost)
363
        This method is called each iteration after popping an index
364
        from the heap, before examining the neighbors.
365

366
        This method can be overloaded to modify the behavior of the MCP
367
        algorithm. An example might be to stop the algorithm when a
368
        certain cumulative cost is reached, or when the front is a
369
        certain distance away from the seed point.
370

371
        This method should return 1 if the algorithm should not check
372
        the current point's neighbors and 2 if the algorithm is now
373
        done.
374
        """
375
        return 0
376

377

378
    cdef void _examine_neighbor(self, INDEX_T index, INDEX_T new_index,
379
                                FLOAT_T offset_length) noexcept:
380
        """ _examine_neighbor(int index, int new_index, float offset_length)
381
        This method is called once for every pair of neighboring nodes,
382
        as soon as both nodes become frozen.
383
        """
384
        pass
385

386

387
    cdef void _update_node(self, INDEX_T index, INDEX_T new_index,
388
                           FLOAT_T offset_length) noexcept:
389
        """ _update_node(int index, int new_index, float offset_length)
390
        This method is called when a node is updated.
391
        """
392
        pass
393

394

395
    def find_costs(self, starts, ends=None, find_all_ends=True,
396
                   max_coverage=1.0, max_cumulative_cost=None, max_cost=None):
397
        """
398
        Find the minimum-cost path from the given starting points.
399

400
        This method finds the minimum-cost path to the specified ending
401
        indices from any one of the specified starting indices. If no end
402
        positions are given, then the minimum-cost path to every position in
403
        the costs array will be found.
404

405
        Parameters
406
        ----------
407
        starts : iterable
408
            A list of n-d starting indices (where n is the dimension of the
409
            `costs` array). The minimum cost path to the closest/cheapest
410
            starting point will be found.
411
        ends : iterable, optional
412
            A list of n-d ending indices.
413
        find_all_ends : bool, optional
414
            If 'True' (default), the minimum-cost-path to every specified
415
            end-position will be found; otherwise the algorithm will stop when
416
            a a path is found to any end-position. (If no `ends` were
417
            specified, then this parameter has no effect.)
418

419
        Returns
420
        -------
421
        cumulative_costs : ndarray
422
            Same shape as the `costs` array; this array records the minimum
423
            cost path from the nearest/cheapest starting index to each index
424
            considered. (If `ends` were specified, not all elements in the
425
            array will necessarily be considered: positions not evaluated will
426
            have a cumulative cost of inf. If `find_all_ends` is 'False', only
427
            one of the specified end-positions will have a finite cumulative
428
            cost.)
429
        traceback : ndarray
430
            Same shape as the `costs` array; this array contains the offset to
431
            any given index from its predecessor index. The offset indices
432
            index into the `offsets` attribute, which is a array of n-d
433
            offsets. In the 2-d case, if offsets[traceback[x, y]] is (-1, -1),
434
            that means that the predecessor of [x, y] in the minimum cost path
435
            to some start position is [x+1, y+1]. Note that if the
436
            offset_index is -1, then the given index was not considered.
437

438
        """
439
        # basic variables to use for end-finding; also fix up the start and end
440
        # lists
441
        cdef BOOL_T use_ends = 0
442
        cdef INDEX_T num_ends = 0
443
        cdef BOOL_T all_ends = find_all_ends
444
        cdef INDEX_T[:] flat_ends
445
        starts = _normalize_indices(starts, self.costs_shape)
446
        if starts is None:
447
            raise ValueError('start points must all be within the costs array')
448
        elif not starts:
449
            raise ValueError('no valid start points to start front' +
450
                             'propagation')
451
        if ends is not None:
452
            ends = _normalize_indices(ends, self.costs_shape)
453
            if ends is None:
454
                raise ValueError('end points must all be within '
455
                                 'the costs array')
456
            use_ends = 1
457
            num_ends = len(ends)
458
            flat_ends = np.array(_ravel_index_fortran(
459
                ends, self.costs_shape), dtype=INDEX_D)
460

461
        # Always perform a reset to (re)initialize our arrays and start
462
        # positions
463
        self._starts, self._ends = starts, ends
464
        self._reset()
465

466
        # Get shorter names for arrays
467
        cdef FLOAT_T[:] flat_costs = self.flat_costs
468
        cdef FLOAT_T[:] flat_cumulative_costs = self.flat_cumulative_costs
469
        cdef OFFSETS_INDEX_T[:] traceback_offsets = self.traceback_offsets
470
        cdef EDGE_T[:, :] flat_pos_edge_map = self.flat_pos_edge_map
471
        cdef EDGE_T[:, :] flat_neg_edge_map = self.flat_neg_edge_map
472
        cdef OFFSET_T[:, :] offsets = self.offsets
473
        cdef INDEX_T[:] flat_offsets = self.flat_offsets
474
        cdef FLOAT_T[:] offset_lengths = self.offset_lengths
475

476
        # Short names for other attributes
477
        cdef heap.FastUpdateBinaryHeap costs_heap = self.costs_heap
478
        cdef DIM_T dim = self.dim
479
        cdef int num_offsets = len(flat_offsets)
480

481
        # Variables used during front propagation
482
        cdef FLOAT_T cost, new_cost, cumcost, new_cumcost, offset_length
483
        cdef INDEX_T index, new_index
484
        cdef BOOL_T is_at_edge, use_offset
485
        cdef INDEX_T d, i, iter
486
        cdef OFFSET_T offset
487
        cdef EDGE_T pos_edge_val, neg_edge_val
488
        cdef int num_ends_found = 0
489
        cdef FLOAT_T inf = np.inf
490
        cdef int goal_reached
491

492
        cdef INDEX_T maxiter = int(max_coverage * flat_costs.size)
493

494
        for iter in range(maxiter):
495

496
            # This is rather like a while loop, except we are guaranteed to
497
            # exit, which is nice during developing to prevent eternal loops.
498

499
            # Find the point with the minimum cost in the heap. Once
500
            # popped, this point's minimum cost path has been found.
501
            if costs_heap.count == 0:
502
                # nothing in the heap: we've found paths to every
503
                # point in the array
504
                break
505

506
            # Get current cumulative cost and index from the heap
507
            cumcost = costs_heap.pop_fast()
508
            index = costs_heap._popped_ref
509

510
            # Record the cost we found to this point
511
            flat_cumulative_costs[index] = cumcost
512

513
            # Check if goal is reached
514
            goal_reached = self.goal_reached(index, cumcost)
515
            if goal_reached > 0:
516
                if goal_reached == 1:
517
                    continue  # Skip neighbors
518
                else:
519
                    break  # Done completely
520

521
            if use_ends:
522
                # If we're only tracing out a path to one or more
523
                # endpoints, check to see if this is an endpoint, and
524
                # if so, if we're done pathfinding.
525
                for i in range(num_ends):
526
                    if index == flat_ends[i]:
527
                        num_ends_found += 1
528
                        break
529
                if (num_ends_found and not all_ends) or \
530
                    num_ends_found == num_ends:
531
                    # if we've found one or all of the end points (as
532
                    # requested), stop searching
533
                    break
534

535
            # Look into the edge map to see if this point is at an
536
            # edge along any axis
537
            is_at_edge = 0
538
            for d in range(dim):
539
                if (flat_pos_edge_map[d, index] != 0 or
540
                    flat_neg_edge_map[d, index] != 0):
541
                    is_at_edge = 1
542
                    break
543

544
            # Now examine the points neighboring the given point
545
            for i in range(num_offsets):
546
                # First, if we're at some edge, scrutinize the offset
547
                # to ensure that it won't put us out-of-bounds. If,
548
                # for example, the edge_map at (x, y) is (-1, 0) --
549
                # though of course we use flat indexing below -- that
550
                # means that (x, y) is along the lower edge of the
551
                # array; thus offsets with -1 or more negative in the
552
                # x-dimension should not be used!
553
                use_offset = 1
554
                if is_at_edge:
555
                    for d in range(dim):
556
                        offset = offsets[i, d]
557
                        pos_edge_val = flat_pos_edge_map[d, index]
558
                        neg_edge_val = flat_neg_edge_map[d, index]
559
                        if (pos_edge_val > 0 and offset >= pos_edge_val) or \
560
                           (neg_edge_val < 0 and offset <= neg_edge_val):
561
                            # the offset puts us out of bounds...
562
                            use_offset = 0
563
                            break
564
                # If not at an edge, or the specific offset doesn't
565
                # push over the edge, then we go on.
566
                if not use_offset:
567
                    continue
568

569
                # using the flat offsets, calculate the new flat index
570
                new_index = index + flat_offsets[i]
571

572
                # Get offset length
573
                offset_length = offset_lengths[i]
574

575
                # If we have already found the best path here then
576
                # ignore this point
577
                if flat_cumulative_costs[new_index] != inf:
578
                    # Give subclass the opportunity to examine these two nodes
579
                    # Note that only when both nodes are "frozen" their
580
                    # cumulative cost is set. By doing the check here, each
581
                    # pair of nodes is checked exactly once.
582
                    self._examine_neighbor(index, new_index, offset_length)
583
                    continue
584

585
                # Get cost and new cost
586
                cost = flat_costs[index]
587
                new_cost = flat_costs[new_index]
588

589
                # If the cost at this point is negative or infinite, ignore it
590
                if new_cost < 0 or new_cost == inf:
591
                    continue
592

593
                # Calculate new cumulative cost
594
                new_cumcost = cumcost + self._travel_cost(cost, new_cost,
595
                                                          offset_length)
596

597
                # Now we ask the heap to append or update the cost to
598
                # this new point, but only if that point isn't already
599
                # in the heap, or it is but the new cost is lower.
600
                # don't push infs into the heap though!
601
                if new_cumcost != inf:
602
                    costs_heap.push_if_lower_fast(new_cumcost, new_index)
603
                    # If we did perform an append or update, we should
604
                    # record the offset from the predecessor to this new
605
                    # point
606
                    if costs_heap._pushed:
607
                        traceback_offsets[new_index] = i
608
                        self._update_node(index, new_index, offset_length)
609

610

611
        # Un-flatten the costs and traceback arrays for human consumption.
612
        cumulative_costs = np.asarray(flat_cumulative_costs)
613
        cumulative_costs = cumulative_costs.reshape(self.costs_shape,
614
                                                    order='F')
615
        traceback = np.asarray(traceback_offsets)
616
        traceback = traceback.reshape(self.costs_shape, order='F')
617
        self.dirty = 1
618
        return cumulative_costs, traceback
619

620

621
    def traceback(self, end):
622
        """traceback(end)
623

624
        Trace a minimum cost path through the pre-calculated traceback array.
625

626
        This convenience function reconstructs the the minimum cost path to a
627
        given end position from one of the starting indices provided to
628
        find_costs(), which must have been called previously. This function
629
        can be called as many times as desired after find_costs() has been
630
        run.
631

632
        Parameters
633
        ----------
634
        end : iterable
635
            An n-d index into the `costs` array.
636

637
        Returns
638
        -------
639
        traceback : list of n-d tuples
640
            A list of indices into the `costs` array, starting with one of
641
            the start positions passed to find_costs(), and ending with the
642
            given `end` index. These indices specify the minimum-cost path
643
            from any given start index to the `end` index. (The total cost
644
            of that path can be read out from the `cumulative_costs` array
645
            returned by find_costs().)
646
        """
647
        if not self.dirty:
648
            raise Exception('find_costs() must be run before traceback()')
649
        ends = _normalize_indices([end], self.costs_shape)
650
        if ends is None:
651
            raise ValueError('the specified end point must be '
652
                             'within the costs array')
653
        traceback = [tuple(ends[0])]
654

655
        cdef INDEX_T flat_position =\
656
             _ravel_index_fortran(ends, self.costs_shape)[0]
657
        if self.flat_cumulative_costs[flat_position] == np.inf:
658
            raise ValueError('no minimum-cost path was found '
659
                             'to the specified end point')
660

661
        # Short names for arrays
662
        cdef OFFSETS_INDEX_T [:] traceback_offsets = self.traceback_offsets
663
        cdef OFFSET_T [:,:] offsets = self.offsets
664
        cdef INDEX_T [:] flat_offsets = self.flat_offsets
665
        # New array
666
        cdef INDEX_T [:] position = np.array(ends[0], dtype=INDEX_D)
667

668
        cdef OFFSETS_INDEX_T offset
669
        cdef DIM_T d
670
        cdef DIM_T dim = self.dim
671
        while 1:
672
            offset = traceback_offsets[flat_position]
673
            if offset == -1:
674
                # At a point where we can go no further: probably a start point
675
                break
676
            flat_position -= flat_offsets[offset]
677
            for d in range(dim):
678
                position[d] -= offsets[offset, d]
679
            traceback.append(tuple(position))
680
        return _reverse(traceback)
681

682

683

684
@cython.boundscheck(False)
685
@cython.wraparound(False)
686
cdef class MCP_Geometric(MCP):
687
    """MCP_Geometric(costs, offsets=None, fully_connected=True)
688

689
    Find distance-weighted minimum cost paths through an n-d costs array.
690

691
    See the documentation for MCP for full details. This class differs from
692
    MCP in that the cost of a path is not simply the sum of the costs along
693
    that path.
694

695
    This class instead assumes that the costs array contains at each position
696
    the "cost" of a unit distance of travel through that position. For
697
    example, a move (in 2-d) from (1, 1) to (1, 2) is assumed to originate in
698
    the center of the pixel (1, 1) and terminate in the center of (1, 2). The
699
    entire move is of distance 1, half through (1, 1) and half through (1, 2);
700
    thus the cost of that move is `(1/2)*costs[1,1] + (1/2)*costs[1,2]`.
701

702
    On the other hand, a move from (1, 1) to (2, 2) is along the diagonal and
703
    is sqrt(2) in length. Half of this move is within the pixel (1, 1) and the
704
    other half in (2, 2), so the cost of this move is calculated as
705
    `(sqrt(2)/2)*costs[1,1] + (sqrt(2)/2)*costs[2,2]`.
706

707
    These calculations don't make a lot of sense with offsets of magnitude
708
    greater than 1. Use the `sampling` argument in order to deal with
709
    anisotropic data.
710
    """
711

712
    def __init__(self, costs, offsets=None, fully_connected=True,
713
                 sampling=None):
714
        """__init__(costs, offsets=None, fully_connected=True, sampling=None)
715

716
        See class documentation.
717
        """
718
        MCP.__init__(self, costs, offsets, fully_connected, sampling)
719
        if np.absolute(self.offsets).max() > 1:
720
            raise ValueError('all offset components must be 0, 1, or -1')
721
        self.use_start_cost = 0
722

723
    cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost,
724
                              FLOAT_T offset_length):
725
        return offset_length * 0.5 * (old_cost + new_cost)
726

727

728

729
@cython.boundscheck(True)
730
@cython.wraparound(True)
731
cdef class MCP_Connect(MCP):
732
    """MCP_Connect(costs, offsets=None, fully_connected=True)
733

734
    Connect source points using the distance-weighted minimum cost function.
735

736
    A front is grown from each seed point simultaneously, while the
737
    origin of the front is tracked as well. When two fronts meet,
738
    create_connection() is called. This method must be overloaded to
739
    deal with the found edges in a way that is appropriate for the
740
    application.
741
    """
742

743
    cdef INDEX_T [:] flat_idmap
744

745

746
    def __init__(self, costs, offsets=None, fully_connected=True,
747
                 sampling=None):
748
        MCP.__init__(self, costs, offsets, fully_connected, sampling)
749

750
        # Create id map to keep track of origin of nodes
751
        self.flat_idmap = np.zeros(self.costs_shape, INDEX_D).ravel('F')
752

753

754
    def _reset(self):
755
        """ Reset the id map.
756
        """
757
        cdef INDEX_T start
758

759
        MCP._reset(self)
760
        starts, _ = self._starts, self._ends
761

762
        # Reset idmap
763
        self.flat_idmap[...] = -1
764
        id = 0
765
        for start in _ravel_index_fortran(starts, self.costs_shape):
766
            self.flat_idmap[start] = id
767
            id += 1
768

769

770
    cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost,
771
                              FLOAT_T offset_length):
772
        """ Equivalent to MCP_Geometric.
773
        """
774
        return offset_length * 0.5 * (old_cost + new_cost)
775

776

777
    cdef void _examine_neighbor(self, INDEX_T index, INDEX_T new_index,
778
                                FLOAT_T offset_length) noexcept:
779
        """ Check whether two fronts are meeting. If so, the flat_traceback
780
        is obtained and a connection is created.
781
        """
782

783
        # Short names
784
        cdef INDEX_T [:] flat_idmap = self.flat_idmap
785
        cdef FLOAT_T [:] flat_cumulative_costs = self.flat_cumulative_costs
786

787
        # Get ids
788
        cdef INDEX_T id1 = flat_idmap[index]
789
        cdef INDEX_T id2 = flat_idmap[new_index]
790

791
        if id2 < 0 or id1 < 0:
792
            pass
793
        elif id2 != id1:
794
            # We reached the 'front' of another seed point!
795
            # Get position/coordinates
796
            pos1, pos2 = _unravel_index_fortran([index, new_index],
797
                                                self.costs_shape)
798
            # Also get the costs, so we can keep the path with the least cost
799
            cost1 = flat_cumulative_costs[index]
800
            cost2 = flat_cumulative_costs[new_index]
801
            # Create connection
802
            self.create_connection(id1, id2, pos1, pos2, cost1, cost2)
803

804

805
    def create_connection(self, id1, id2, tb1, tb2, cost1, cost2):
806
        """ create_connection id1, id2, pos1, pos2, cost1, cost2)
807

808
        Overload this method to keep track of the connections that are
809
        found during MCP processing. Note that a connection with the
810
        same ids can be found multiple times (but with different
811
        positions and costs).
812

813
        At the time that this method is called, both points are "frozen"
814
        and will not be visited again by the MCP algorithm.
815

816
        Parameters
817
        ----------
818
        id1 : int
819
            The seed point id where the first neighbor originated from.
820
        id2 : int
821
            The seed point id where the second neighbor originated from.
822
        pos1 : tuple
823
            The index of of the first neighbor in the connection.
824
        pos2 : tuple
825
            The index of of the second neighbor in the connection.
826
        cost1 : float
827
            The cumulative cost at `pos1`.
828
        cost2 : float
829
            The cumulative costs at `pos2`.
830
        """
831
        pass
832

833

834
    cdef void _update_node(self, INDEX_T index, INDEX_T new_index,
835
                           FLOAT_T offset_length) noexcept:
836
        """ Keep track of the id map so that we know which seed point
837
        a certain front originates from.
838
        """
839
        self.flat_idmap[new_index] = self.flat_idmap[index]
840

841

842

843
@cython.boundscheck(False)
844
@cython.wraparound(False)
845
cdef class MCP_Flexible(MCP):
846
    """MCP_Flexible(costs, offsets=None, fully_connected=True)
847

848
    Find minimum cost paths through an N-d costs array.
849

850
    See the documentation for MCP for full details. This class differs from
851
    MCP in that several methods can be overloaded (from pure Python) to
852
    modify the behavior of the algorithm and/or create custom algorithms
853
    based on MCP. Note that goal_reached can also be overloaded in the
854
    MCP class.
855

856
    """
857

858
    def travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost,
859
                    FLOAT_T offset_length):
860
        """ travel_cost(old_cost, new_cost, offset_length)
861
        This method calculates the travel cost for going from the
862
        current node to the next. The default implementation returns
863
        new_cost. Overload this method to adapt the behaviour of the
864
        algorithm.
865
        """
866
        return new_cost
867

868

869
    def examine_neighbor(self, INDEX_T index, INDEX_T new_index,
870
                         FLOAT_T offset_length):
871
        """ examine_neighbor(index, new_index, offset_length)
872
        This method is called once for every pair of neighboring nodes,
873
        as soon as both nodes are frozen.
874

875
        This method can be overloaded to obtain information about
876
        neighboring nodes, and/or to modify the behavior of the MCP
877
        algorithm. One example is the MCP_Connect class, which checks
878
        for meeting fronts using this hook.
879
        """
880
        pass
881

882

883
    def update_node(self, INDEX_T index, INDEX_T new_index,
884
                    FLOAT_T offset_length):
885
        """ update_node(index, new_index, offset_length)
886
        This method is called when a node is updated, right after
887
        new_index is pushed onto the heap and the traceback map is
888
        updated.
889

890
        This method can be overloaded to keep track of other arrays
891
        that are used by a specific implementation of the algorithm.
892
        For instance the MCP_Connect class uses it to update an id map.
893
        """
894
        pass
895

896

897
    cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost,
898
                              FLOAT_T offset_length):
899
        return self.travel_cost(old_cost, new_cost, offset_length)
900

901

902
    cdef void _examine_neighbor(self, INDEX_T index, INDEX_T new_index,
903
                                FLOAT_T offset_length) noexcept:
904
        self.examine_neighbor(index, new_index, offset_length)
905

906

907
    cdef void _update_node(self, INDEX_T index, INDEX_T new_index,
908
                           FLOAT_T offset_length) noexcept:
909
        self.update_node(index, new_index, offset_length)
910