~wecacuee/mrol/mrol-dev

« back to all changes in this revision

Viewing changes to mrol_mapping/bresenhamline.py

  • Committer: Vikas Dhiman
  • Date: 2012-05-22 15:26:38 UTC
  • Revision ID: wecacuee@gmail.com-20120522152638-j6srlkgrpdq4kdfe
Added support for userdata in terms of colored voxels. Added benchmark test and corresponding data

Show diffs side-by-side

added added

removed removed

Lines of Context:
2
2
N-D Bresenham line algo
3
3
"""
4
4
import numpy as np
5
 
def bresenhamline_nslope(slope):
 
5
def _bresenhamline_nslope(slope):
6
6
    """
7
7
    Normalize slope for Bresenham's line algorithm.
8
8
 
9
9
    >>> s = np.array([[-2, -2, -2, 0]])
10
 
    >>> bresenhamline_nslope(s)
 
10
    >>> _bresenhamline_nslope(s)
11
11
    array([[-1., -1., -1.,  0.]])
12
12
 
13
13
    >>> s = np.array([[0, 0, 0, 0]])
14
 
    >>> bresenhamline_nslope(s)
 
14
    >>> _bresenhamline_nslope(s)
15
15
    array([[ 0.,  0.,  0.,  0.]])
16
16
 
17
17
    >>> s = np.array([[0, 0, 9, 0]])
18
 
    >>> bresenhamline_nslope(s)
 
18
    >>> _bresenhamline_nslope(s)
19
19
    array([[ 0.,  0.,  1.,  0.]])
20
20
    """
21
21
    scale = np.amax(np.abs(slope), axis=1).reshape(-1, 1)
22
22
    zeroslope = (scale == 0).all(1)
23
23
    scale[zeroslope] = np.ones(1)
24
 
    normalizedslope = np.array(slope, dtype=np.double) / scale
 
24
    normalizedslope = slope.astype(np.double) / scale
25
25
    normalizedslope[zeroslope] = np.zeros(slope[0].shape)
26
26
    return normalizedslope
27
27
 
28
 
def bresenhamline_step(start, normalizedslope, err):
29
 
    """
30
 
    Performs a single incremental step according to Bresenham's algorithm.
31
 
 
32
 
    >>> s = np.array([[4, 4, 4, 0]])
33
 
    >>> ns = bresenhamline_nslope(np.array([[-2, -2, -2, 0]]))
34
 
    >>> e = np.zeros((1, 4))
35
 
    >>> bresenhamline_step(s, ns, e)
36
 
    (array([[ 3.,  3.,  3.,  0.]]), array([[ 0.,  0.,  0.,  0.]]))
37
 
 
38
 
    >>> s = np.array([[0, 0, 9, 0]])
39
 
    >>> ns = bresenhamline_nslope(-1 * s)
40
 
    >>> e = np.zeros((1, 4))
41
 
    """
42
 
    nextpoint = start + normalizedslope + err
43
 
    nextvox = np.rint(nextpoint)
44
 
    err = nextpoint - nextvox
45
 
    return nextvox, err
46
 
 
47
 
def bresenhamline_gen(start, end):
48
 
    """
49
 
    Returns a generator of points from (start, end] by ray tracing a line b/w
50
 
    the points.
51
 
    Parameters:
52
 
        start: An array of start points (number of points x dimension)
53
 
        end:   An end points (1 x dimension)
54
 
            or An array of end point corresponding to each start point
55
 
                (number of points x dimension)
56
 
 
57
 
    Returns:
58
 
        A generator that returns a point being traversed at each step
59
 
        corresponding to each start point.
60
 
 
61
 
    >>> [p for p in bresenhamline_gen(np.array([[3, 1, 3, 0]]), np.zeros(4))]
62
 
    [array([[2, 1, 2, 0]]), array([[1, 0, 1, 0]]), array([[0, 0, 0, 0]])]
63
 
    """
64
 
    dim = start.shape[1]
65
 
    err = np.zeros(start.shape)
66
 
    slope = np.array(end - start, dtype=np.float32)
67
 
    nslope = bresenhamline_nslope(slope)
68
 
    cur_vox = start
69
 
    while (np.sum(nslope != 0)):
70
 
        cur_vox, err = bresenhamline_step(cur_vox, nslope, err)
71
 
        # reached end ?
72
 
        nslope[(np.abs(end - cur_vox) < 1).all(1)] = np.zeros(dim) 
73
 
        yield np.array(cur_vox, dtype=start.dtype)
 
28
def _bresenhamlines(start, end, max_iter):
 
29
    """
 
30
    Returns npts lines of length max_iter each. (npts x max_iter x dimension) 
 
31
 
 
32
    >>> s = np.array([[3, 1, 9, 0],[0, 0, 3, 0]])
 
33
    >>> _bresenhamlines(s, np.zeros(s.shape[1]), max_iter=-1)
 
34
    array([[[ 3,  1,  8,  0],
 
35
            [ 2,  1,  7,  0],
 
36
            [ 2,  1,  6,  0],
 
37
            [ 2,  1,  5,  0],
 
38
            [ 1,  0,  4,  0],
 
39
            [ 1,  0,  3,  0],
 
40
            [ 1,  0,  2,  0],
 
41
            [ 0,  0,  1,  0],
 
42
            [ 0,  0,  0,  0]],
 
43
    <BLANKLINE>
 
44
           [[ 0,  0,  2,  0],
 
45
            [ 0,  0,  1,  0],
 
46
            [ 0,  0,  0,  0],
 
47
            [ 0,  0, -1,  0],
 
48
            [ 0,  0, -2,  0],
 
49
            [ 0,  0, -3,  0],
 
50
            [ 0,  0, -4,  0],
 
51
            [ 0,  0, -5,  0],
 
52
            [ 0,  0, -6,  0]]])
 
53
    """
 
54
    if max_iter == -1:
 
55
        max_iter = np.amax(np.amax(np.abs(end - start), axis=1))
 
56
    npts, dim = start.shape
 
57
    nslope = _bresenhamline_nslope(end - start)
 
58
 
 
59
    # steps to iterate on
 
60
    stepseq = np.arange(1, max_iter + 1)
 
61
    stepmat = np.tile(stepseq, (dim, 1)).T
 
62
 
 
63
    # some hacks for broadcasting properly
 
64
    bline = start[:, np.newaxis, :] + nslope[:, np.newaxis, :] * stepmat
 
65
 
 
66
    # Approximate to nearest int
 
67
    return np.rint(bline).astype(start.dtype)
74
68
 
75
69
def bresenhamline(start, end, max_iter=5):
76
70
    """
81
75
        end:   An end points (1 x dimension)
82
76
            or An array of end point corresponding to each start point
83
77
                (number of points x dimension)
84
 
        max_iter: Max points to traverse
 
78
        max_iter: Max points to traverse. if -1, maximum number of required
 
79
                  points are traversed
85
80
 
86
81
    Returns:
87
82
        linevox (n x dimension) A cumulative array of all points traversed by
88
83
        all the lines so far.
89
84
 
90
 
    >>> s = np.array([[0, 0, 3, 0]])
91
 
    >>> bresenhamline(s, np.zeros(s.shape[1]), max_iter=9)
92
 
    array([[0, 0, 2, 0],
93
 
           [0, 0, 1, 0],
94
 
           [0, 0, 0, 0]])
95
 
 
96
 
    >>> s = np.array([[3, 1, 9, 0]])
97
 
    >>> bresenhamline(s, np.zeros(s.shape[1]), max_iter=9)
98
 
    array([[3, 1, 8, 0],
99
 
           [2, 1, 7, 0],
100
 
           [2, 1, 6, 0],
101
 
           [2, 1, 5, 0],
102
 
           [1, 0, 4, 0],
103
 
           [1, 0, 3, 0],
104
 
           [1, 0, 2, 0],
105
 
           [0, 0, 1, 0],
106
 
           [0, 0, 0, 0]])
 
85
    >>> s = np.array([[3, 1, 9, 0],[0, 0, 3, 0]])
 
86
    >>> bresenhamline(s, np.zeros(s.shape[1]), max_iter=-1)
 
87
    array([[ 3,  1,  8,  0],
 
88
           [ 2,  1,  7,  0],
 
89
           [ 2,  1,  6,  0],
 
90
           [ 2,  1,  5,  0],
 
91
           [ 1,  0,  4,  0],
 
92
           [ 1,  0,  3,  0],
 
93
           [ 1,  0,  2,  0],
 
94
           [ 0,  0,  1,  0],
 
95
           [ 0,  0,  0,  0],
 
96
           [ 0,  0,  2,  0],
 
97
           [ 0,  0,  1,  0],
 
98
           [ 0,  0,  0,  0],
 
99
           [ 0,  0, -1,  0],
 
100
           [ 0,  0, -2,  0],
 
101
           [ 0,  0, -3,  0],
 
102
           [ 0,  0, -4,  0],
 
103
           [ 0,  0, -5,  0],
 
104
           [ 0,  0, -6,  0]])
107
105
    """
108
 
    dim = start.shape[1]
109
 
    linevox = np.zeros((0, dim), dtype=start.dtype)
110
 
    linegen = bresenhamline_gen(start, end)
111
 
    try:
112
 
        for i in range(max_iter):
113
 
            cur_vox = linegen.next()
114
 
            linevox = np.vstack((linevox, cur_vox))
115
 
    except StopIteration:
116
 
        pass
117
 
    return np.array(linevox, dtype=start.dtype)
 
106
    # Return the points as a single array
 
107
    return _bresenhamlines(start, end, max_iter).reshape(-1, start.shape[-1])
 
108
 
118
109
 
119
110
if __name__ == "__main__":
120
111
    import doctest