~wecacuee/mrol/mrol-dev

« back to all changes in this revision

Viewing changes to mrol_mapping/poseutil.py

  • Committer: Julian Ryde
  • Date: 2011-06-03 16:52:15 UTC
  • Revision ID: julian_ryde-20110603165215-oqr6yo9nwng3j9nb
Initial import from revision 206 of Julian Ryde's mrol_mapping repository.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
import numpy as np
 
2
import quantizer
 
3
 
 
4
# TODO replace a transformations  library?
 
5
# transformations.py now in lib directory
 
6
# This library also has random rotation pose generation 
 
7
 
 
8
# TODO make all interfaces consistent and using Pose3D objects
 
9
 
 
10
# TODO need to test mat and posetuple functions properly rather than just testing that 
 
11
# they invert one another
 
12
 
 
13
def mat(pose):
 
14
    '''Convert pose tuple in form x, y, z, rotx, roty, rotz to matrix'''
 
15
    #M = transformations.euler_matrix(pose[3], pose[4], pose[5], 'sxyz')
 
16
    #M[:3, 3] = pose[:3]
 
17
    P = Pose3D(pose)
 
18
    return P.getMatrix()
 
19
 
 
20
def posetuple(posemat):
 
21
    '''Convert pose matrix 4x4 to pose tuple x, y, z, rotx, roty, rotz'''
 
22
    P = Pose3D()
 
23
    P.setMatrix(posemat)
 
24
    return P.get()
 
25
    
 
26
# TODO combine tuplepose and posetuple into one function
 
27
def tuplepose(pose):
 
28
    '''Convert an array like 6 element pose representation to a tuple x, y, z, rotx, roty, rotz'''
 
29
    return (pose[0], pose[1], pose[2], pose[3], pose[4], pose[5])
 
30
    
 
31
# TODO enable transformPoints function to take either matrix or pose tuple, and modify 
 
32
# where this function is invoked accordingly
 
33
def transformPoints(xyz, posemat):
 
34
    '''returns the transform of a list of points n by 3 by applying the 
 
35
    matrix transform described by mat'''
 
36
    # TODO probably would be faster to transform in place
 
37
    #M = np.vstack((xyz[:, 0:3].T, np.ones((1, xyz.shape[0]))))
 
38
    #return np.dot(mat, M).T[:, 0:3]
 
39
    # below is faster?
 
40
    rotated = np.dot(posemat[:3, :3], xyz.T).T
 
41
    rotated += posemat[:3, 3] 
 
42
    return rotated
 
43
    
 
44
# TODO combine transformPointsTuple into transformPoints
 
45
def transformPointsTuple(xyz, pose_tuple):
 
46
    pose_mat = mat(pose_tuple)
 
47
    return transformPoints(xyz, pose_mat)
 
48
 
 
49
def refineposes(poses, deltaposition, deltaorientation):
 
50
    '''Produces a new set of poses that are around the supplied poses'''
 
51
    # TODO this needs to be done properly for 3D
 
52
    # TODO does not handle a single pose
 
53
   
 
54
    # each pose should spawn 3 new poses for each dimension
 
55
    # this definitely produces better results in 2D tested by testing on 2D mapping 
 
56
    # and looking at the voxel count of the final map
 
57
    #mg = np.mgrid[-1:2, -1:2, -1:2, -1:2, -1:2, -1:2].astype(float)
 
58
    # each pose should spawn 2 new poses for each dimension
 
59
 
 
60
    mg = np.mgrid[-1:1, -1:1, -1:1, -1:1, -1:1, -1:1].astype(float)
 
61
    mg += 0.5
 
62
 
 
63
    mg.shape = 6, -1
 
64
    mg = mg.T
 
65
 
 
66
    # TODO remove this nasty hack for 3DOF speed
 
67
    #mg[:, [2, 3, 4]] = 0
 
68
    #mg = quantizer.uniquerows(mg, 3)
 
69
 
 
70
    mg[:, 0:3] *= deltaposition
 
71
    mg[:, 3:6] *= deltaorientation
 
72
    #zerorows = np.zeros((6, mg.shape[1]))
 
73
    #mg = np.vstack((mg[0:3, :]*deltaposition, zerorows, mg[3:6, :]*deltaorientation))
 
74
 
 
75
    refinedposes = np.empty((poses.shape[0]*mg.shape[0], 6))
 
76
    row = 0
 
77
    for p in poses:
 
78
        ptile = np.tile(p, (mg.shape[0], 1))
 
79
        newposes = ptile + mg
 
80
        refinedposes[row:row+mg.shape[0], :] = newposes
 
81
        row += mg.shape[0]
 
82
 
 
83
    # TODO remove those outside of centrepose + spreadpose
 
84
    return quantizer.uniquerows(refinedposes, 3)
 
85
 
 
86
def uniformposes(centrepose, spreadpose, positionres, orientationres):
 
87
    """ Produce uniform distribution of poses with positions separated by 
 
88
    positionres and orientations separated by orientationres, centered around 
 
89
    center pose and spanning +/- spreadpose """
 
90
    # TODO this needs to be done properly for 3D poses
 
91
    # see papers on generating uniform random 3D rotations
 
92
    cp = np.array(centrepose.get())
 
93
    sp = np.array(spreadpose.get())
 
94
    B = cp - sp 
 
95
    E = cp + sp + 0.0001 
 
96
    # this small addition is to fix when the spreadpose is 0 still want to 
 
97
    # generate some poses rather than none
 
98
    S = np.hstack((np.tile(positionres, 3), np.tile(orientationres, 3)))
 
99
    pose_count = np.prod((E-B)/S)
 
100
    assert pose_count < 1e5, 'Requested too many poses %d' % int(pose_count)
 
101
 
 
102
    mg = np.mgrid[B[0]:E[0]:S[0], B[1]:E[1]:S[1], B[2]:E[2]:S[2], B[3]:E[3]:S[3], B[4]:E[4]:S[4], B[5]:E[5]:S[5]]
 
103
    mg.shape = 6, -1
 
104
    return quantizer.uniquerows(mg.T, 3)
 
105
    
 
106
def randomrot():
 
107
    """ Return an uniformly distributed random rotation matrix """
 
108
    # TODO write a test for this somehow. NFI.
 
109
    # make a random quaternion 
 
110
    randelements = np.random.rand(3)
 
111
    lin1 = np.sqrt(1.0 - randelements[0])
 
112
    lin2 = np.sqrt(randelements[0])
 
113
    rot1 = 2*np.pi * randelements[1]
 
114
    rot2 = 2*np.pi * randelements[2]
 
115
    x = np.sin(rot1)*lin1
 
116
    y = np.cos(rot1)*lin1
 
117
    z = np.sin(rot2)*lin2
 
118
    w = np.cos(rot2)*lin2
 
119
    quat = np.array([x,y,z,w])
 
120
    quat *= np.sqrt(2.0/np.dot(quat,quat)) # not sure if we need this line
 
121
    q = np.outer(quat,quat) # assumes quat is properly normalised... hrmm
 
122
    mat = np.array([[1.0-q[1, 1]-q[2, 2],     q[0, 1]-q[2, 3],              q[0, 2]+q[1, 3]], 
 
123
                    [    q[0, 1]+q[2, 3],     1.0-q[0, 0]-q[2, 2],          q[1, 2]-q[0, 3]], 
 
124
                    [    q[0, 2]-q[1, 3],         q[1, 2]+q[0, 3],      1.0-q[0, 0]-q[1, 1]]])
 
125
    return mat
 
126
    
 
127
    
 
128
 
 
129
# TODO allow pose_str function to also take in a pose in matrix form
 
130
def pose_str(P):
 
131
    '''Convert 6 element pose x, y, z, rotx, roty, rotz to a string with angles 
 
132
    in degrees'''
 
133
    pose = Pose3D(P)
 
134
    return str(pose)
 
135
    
 
136
def check_and_transform_points(xyzs,pose):
 
137
    '''Check the type of pose passed and perform the transformation'''
 
138
    if type(pose) == type(Pose3D()):
 
139
        if pose == Pose3D():
 
140
            xyzs = xyzs # do nothing
 
141
        else:
 
142
            xyzs = transformPoints(xyzs, pose.getMatrix())   
 
143
        pose_return = pose
 
144
    else:
 
145
        if pose == None or (np.array(pose) == np.array((0, 0, 0, 0, 0, 0))).all() or (np.array(pose) == np.array((0.0, 0.0, 0.0, 0.0, 0.0, 0.0))).all():
 
146
            xyzs = xyzs # do nothing
 
147
            pose_return = Pose3D()
 
148
        else:
 
149
            xyzs = transformPointsTuple(xyzs, pose)
 
150
            pose_return = Pose3D(pose)
 
151
    return xyzs, pose_return
 
152
    
 
153
def inverse(pose):
 
154
    '''compute the inverse transform of that given'''
 
155
    posemat = mat(pose)
 
156
    invPose = np.matrix(posemat).I
 
157
    return posetuple(invPose)
 
158
    
 
159
# TODO use an array to store pose rather than self.x, self.y, ...
 
160
class Pose3D:
 
161
    """ 3D pose """
 
162
    def __init__(self, X=(0, 0, 0, 0, 0, 0)):
 
163
        """ Initialises a pose with no translation and no rotation """
 
164
        #if type(posestr) == str:
 
165
        #    fields = posestr.split()
 
166
        #else:
 
167
        #    fields = posestr
 
168
        
 
169
        # TODO Something like below so we can assign from one Pose3D to another
 
170
        #if type(X) == type(Pose3D()):
 
171
        #    X = X.getTuple()
 
172
        
 
173
        self.set(X)
 
174
        self.tol = 6 # number of decimals places to round to
 
175
 
 
176
    def __repr__(self):
 
177
        return self.__str__()
 
178
 
 
179
    def __str__(self):
 
180
        return '%.4f %.4f %.4f %.4f %.4f %.4f' % (self.x, self.y, self.z, np.degrees(self.rotx), np.degrees(self.roty), np.degrees(self.rotz))
 
181
 
 
182
    # __eq__ and __hash__ must be consistent
 
183
    # TODO fix to handle equal tolerance maybe with rounding
 
184
    def __eq__(self, o):
 
185
        if o == None:
 
186
            return False
 
187
        return np.allclose(self.get(), o.get(), atol=0.0001)
 
188
    
 
189
    def __hash__(self):
 
190
        X = np.around(self.get(), self.tol)
 
191
        return hash(tuple(X))
 
192
 
 
193
    def is_close(self, o, delta_pos, delta_angle_rad):
 
194
        '''Returns true if the pose is within delta_pos position and 
 
195
        delta_angle orientation, delta_angle is in radians'''
 
196
        # TODO will be better when the pose components are stored as a numpy 
 
197
        # array internally 
 
198
        if abs(self.x - o.x) > delta_pos: return False
 
199
        if abs(self.y - o.y) > delta_pos: return False
 
200
        if abs(self.z - o.z) > delta_pos: return False
 
201
        if abs(self.rotx - o.rotx) > delta_angle_rad: return False
 
202
        if abs(self.roty - o.roty) > delta_angle_rad: return False
 
203
        if abs(self.rotz - o.rotz) > delta_angle_rad: return False
 
204
        return True
 
205
 
 
206
    def copy(self):
 
207
        c = Pose3D()
 
208
        c.setMatrix(self.getMatrix())
 
209
        return c
 
210
 
 
211
    def inverse(self):
 
212
        """ returns a pose which is the inverse of this pose 
 
213
        """
 
214
        M = np.linalg.inv(self.getMatrix())
 
215
        p = Pose3D()
 
216
        p.setMatrix(M)
 
217
        return p
 
218
 
 
219
    def set_position(self, position):
 
220
        self.x = position[0]
 
221
        self.y = position[1]
 
222
        self.z = position[2]
 
223
 
 
224
    def set(self, X):
 
225
        X = np.ravel(X) # handle lists, tuples and 6x1 and 1x6 arrays
 
226
        self.x = float(X[0])
 
227
        self.y = float(X[1])
 
228
        self.z = float(X[2])
 
229
        self.rotx = float(X[3])
 
230
        self.roty = float(X[4])
 
231
        self.rotz = float(X[5])
 
232
 
 
233
    def get(self):
 
234
        return (self.x, self.y, self.z, self.rotx, self.roty, self.rotz)
 
235
 
 
236
    # TODO refactor these identical functions together?
 
237
    def getTuple(self):
 
238
        return (self.x, self.y, self.z, self.rotx, self.roty, self.rotz)
 
239
 
 
240
    def getMatrix(self):
 
241
        """ returns 4x4 matrix representing this pose """
 
242
        sq = np.sin(self.rotx)
 
243
        cq = np.cos(self.rotx)
 
244
        Rx = [[1, 0, 0], [0, cq, -sq], [0, sq, cq]]
 
245
        sq = np.sin(self.roty)
 
246
        cq = np.cos(self.roty)
 
247
        Ry = [[cq, 0, sq], [0, 1, 0], [-sq, 0, cq]]
 
248
        sq = np.sin(self.rotz)
 
249
        cq = np.cos(self.rotz)
 
250
        Rz = [[cq, -sq, 0], [sq, cq, 0], [0, 0, 1]]
 
251
        # this order is important
 
252
        #R = np.dot(Rx, np.dot(Ry, Rz))
 
253
        R = np.dot(Rz, np.dot(Ry, Rx))
 
254
        T = [[self.x], [self.y], [self.z]]
 
255
        trans = np.hstack((R, T))
 
256
        trans = np.vstack([trans, [0, 0, 0, 1]])
 
257
        return trans.astype(np.single)
 
258
 
 
259
    def setMatrix(self, mat):
 
260
        """ sets the pose according to the 4x4 transform matrix """
 
261
        self.rotx = np.arctan2(mat[2, 1], mat[2, 2]);
 
262
        self.roty = np.arcsin(-mat[2, 0]);
 
263
        self.rotz = np.arctan2(mat[1, 0], mat[0, 0]);
 
264
        self.x = mat[0, 3]
 
265
        self.y = mat[1, 3]
 
266
        self.z = mat[2, 3]
 
267
 
 
268
    def getCombined(self, p):
 
269
        A = self.getMatrix()
 
270
        B = p.getMatrix()
 
271
        C = Pose3D()
 
272
        C.setMatrix(np.dot(A, B))
 
273
        return C
 
274
 
 
275
    def transformPoints(self, xyz):
 
276
        """ returns a transform of a list of points n by 3 by applying the 
 
277
        matrix transform described by this pose """
 
278
        return transformPoints(xyz, self.getMatrix())
 
279
        
 
280
    def inverse(self):
 
281
        return Pose3D(inverse(self.getTuple()))
 
282
        
 
283
def axis_order(horz_rot, vert_rot, axis):
 
284
    assert axis in [0,1,2], "axis must be 0, 1 or 2 (x, y or z)"
 
285
    if axis == 2:
 
286
        return ([horz_rot, vert_rot, 0])
 
287
    elif axis == 1:
 
288
        return ([horz_rot, 0, vert_rot])
 
289
    elif axis == 0:    
 
290
        return ([0, horz_rot, vert_rot])
 
291
        
 
292
def scan_pose_list(FOVV, resV, FOVH, resH, axis=0):
 
293
    ''' Create a list of poses that cover the angular region defined by 
 
294
    the input arguments, about the defined axis'''
 
295
    assert axis in [0,1,2], "axis must be 0, 1 or 2 (x, y or z)"
 
296
    print "scan poses based on:"
 
297
    print FOVV, resV, FOVH, resH
 
298
    assert np.remainder(np.round(FOVV/resV),2) == 0, "FOVV/resV must be an even integer"
 
299
    assert np.remainder(np.round(FOVH/resH),2) == 0, "FOVH/resH must be an even integer"
 
300
    numV = 2*np.ceil(FOVV/resV/2)+1
 
301
    numH = 2*np.ceil(FOVH/resH/2)+1
 
302
    # must be odd so we have a definitive center pixel
 
303
    if np.remainder(numV,2) != 1:
 
304
        numV = numV + 1
 
305
    if np.remainder(numV,2) != 1:
 
306
        numV = numV + 1
 
307
    numposes = numV*numH
 
308
    center = np.floor(numposes/2)
 
309
    poses_ang = np.zeros((numposes,3))
 
310
    poses_lin = np.zeros((numposes,3))
 
311
    # top left
 
312
    count = 0
 
313
    for i in np.arange(-FOVH/2,0,resH):
 
314
        horz_rot = i
 
315
        for j in np.arange(0,FOVV/2+resV,resV):
 
316
            vert_rot = j
 
317
            poses_ang[count,:] = axis_order(horz_rot, vert_rot, axis)
 
318
            count = count + 1
 
319
        vert_rot = 0
 
320
    # top right, center pixel
 
321
    for i in np.arange(0,FOVH/2+resH,resH):
 
322
        horz_rot = i
 
323
        for j in np.arange(0,FOVV/2+resV,resV):
 
324
            vert_rot = j
 
325
            poses_ang[count,:] = axis_order(horz_rot, vert_rot, axis)
 
326
            count = count + 1
 
327
        vert_rot = 0
 
328
    # bottom left
 
329
    for i in np.arange(-FOVH/2,0,resH):
 
330
        horz_rot = i
 
331
        for j in np.arange(-FOVV/2,0,resV):
 
332
            vert_rot = j
 
333
            poses_ang[count,:] = axis_order(horz_rot, vert_rot, axis)
 
334
            count = count + 1
 
335
        vert_rot = 0
 
336
    # bottom right
 
337
    for i in np.arange(0,FOVH/2+resH,resH):
 
338
        horz_rot = i
 
339
        for j in np.arange(-FOVV/2,0,resV):
 
340
            vert_rot = j
 
341
            poses_ang[count,:] = axis_order(horz_rot, vert_rot, axis)
 
342
            count = count + 1
 
343
        vert_rot = 0
 
344
    poses = np.hstack([poses_lin,poses_ang])
 
345
    return poses