dwm1000: Always send transform even if there are no valid sensors values
[ros_wild_thumper.git] / scripts / dwm1000.py
1 #!/usr/bin/env python
2 # -*- coding: iso-8859-15 -*-
3
4 VISULAIZE = False
5
6 import threading
7 import struct
8 import rospy
9 import tf
10 import numpy as np
11 from math import *
12 from datetime import datetime
13 from i2c import i2c
14 from time import sleep
15 from std_msgs.msg import Float32
16 from nav_msgs.msg import Odometry
17 from wild_thumper.srv import DWM1000Center, DWM1000CenterResponse
18 if VISULAIZE:
19         import matplotlib.pyplot as plt
20
21 class simple_kalman:
22         def __init__(self, x_est, P_est, Q, R):
23                 self.x_est = x_est # Systemzustand
24                 self.P_est = P_est # Fehlerkovarianz
25                 self.Q = Q # Systemrauschen
26                 self.R = R # Varianz des Messfehlers
27
28         def run(self, y):
29                 # Korrektur mit der Messung
30                 # (1) Berechnung der Kalman Verstärkung
31                 K = self.P_est/(self.R + self.P_est)
32                 # (2) Korrektur der Schätzung mit der Messung y
33                 x = self.x_est + K*(y - self.x_est)
34                 # (3) Korrektur der Fehlerkovarianzmatrix
35                 P = (1-K)*self.P_est
36                 #
37                 # Prädiktion
38                 # (1) Prädiziere den Systemzustand
39                 self.x_est = x
40                 # (2) Präzidiere die Fehlerkovarianzmatrix
41                 self.P_est = P + self.Q
42
43                 return x
44
45         def set_measure_cov(self, R):
46                 if R > 0:
47                         self.R = R
48
49 class DW1000(threading.Thread):
50         def __init__(self, name, addr, offset):
51                 threading.Thread.__init__(self)
52                 self.setDaemon(1)
53                 self.dist = 0
54                 self.offset = offset
55                 self.addr = addr
56                 self.name = name
57                 self.last_update = datetime.min
58
59                 self.pub = rospy.Publisher(name, Float32, queue_size=16)
60
61                 self.start()
62
63         def get_value(self):
64                 dev = i2c(self.addr)
65                 ret = struct.unpack("f", dev.read(4))
66                 dev.close()
67                 return ret[0]
68
69         def distance(self):
70                 return self.dist
71
72         # Returns each distance only if current
73         def distance_valid(self):
74                 if (datetime.now() - self.last_update).seconds < 1:
75                         return self.dist
76                 return None
77
78         def run(self):
79                 last_val = 10
80                 while not rospy.is_shutdown():
81                         val = self.get_value()
82                         if abs(val - last_val)  > 50:
83                                 rospy.logwarn("Ignoring values too far apart %s: %.2f - %.2f", self.name, val, last_val)
84                         elif not isnan(val):
85                                 self.dist = val + self.offset
86                                 self.last_update = datetime.now()
87                                 self.pub.publish(self.distance())
88                                 last_val = val
89                         sleep(0.1)
90
91 class Position:
92         def __init__(self):
93                 # Varianz des Messfehlers
94                 Rx = 0.2
95                 Ry = 0.05
96                 # Fehlerkovarianz
97                 P_est_x = 0.02
98                 P_est_y = 0.01
99                 # Systemrauschen
100                 Qx = 0.001
101                 Qy = 0.00001
102                 # TODO: Replace 2x 1d kalman with 2d kalman
103                 # so the polyfit variance can be replaced with a covariance matrix
104                 self.filter_x = simple_kalman(1.0, P_est_x, Qx, Rx)
105                 self.filter_y = simple_kalman(0.0, P_est_y, Qy, Ry)
106                 self.speed_x = 0
107                 self.speed_y = 0
108                 self.speed_z = 0
109                 self.last_time = rospy.Time.now()
110                 rospy.Subscriber("/odom_combined", Odometry, self.odomReceived)
111
112         def odomReceived(self, msg):
113                 self.speed_x = msg.twist.twist.linear.x
114                 self.speed_y = msg.twist.twist.linear.y
115                 self.speed_z = msg.twist.twist.angular.z
116
117         def filter(self, x, y):
118                 # Correct estimation with speed
119                 current_time = rospy.Time.now()
120                 dt = (current_time - self.last_time).to_sec()
121                 # Subtract vehicle speed
122                 pos = np.array([self.filter_x.x_est, self.filter_y.x_est])
123                 # translation
124                 pos -= np.array([self.speed_x*dt, self.speed_y*dt])
125                 # rotation
126                 rot = np.array([[np.cos(self.speed_z*dt), -np.sin(self.speed_z*dt)],
127                                 [np.sin(self.speed_z*dt),  np.cos(self.speed_z*dt)]])
128                 pos = np.dot(pos, rot)
129                 # copy back
130                 self.filter_x.x_est = pos[0]
131                 self.filter_y.x_est = pos[1]
132
133                 # run kalman if new measurements are valid
134                 if x != None and y != None:
135                         x = self.filter_x.run(x)
136                         y = self.filter_y.run(y)
137
138                         # Update variance
139                         dist = np.linalg.norm([x, y])
140                         self.filter_x.set_measure_cov(np.polyval([0.018, 0.0, 0.0], dist))
141                         self.filter_y.set_measure_cov(np.polyval([0.006, 0.0, 0.0], dist))
142                 else:
143                         x = self.filter_x.x_est
144                         y = self.filter_y.x_est
145
146                 self.last_time = current_time
147                 return x,y
148
149
150 def handle_center_call(req):
151         diff = dwleft.distance_valid() - dwright.distance_valid()
152         dwleft.offset -= diff/2
153         dwright.offset += diff/2
154         rospy.loginfo("Centering to %.2f %.2f", dwleft.offset, dwright.offset)
155         return DWM1000CenterResponse()
156
157 if __name__ == "__main__":
158         rospy.init_node('DWM1000', log_level=rospy.DEBUG)
159         dwleft  = DW1000("uwb_dist_left",  0xc2, +0.0)
160         dwright = DW1000("uwb_dist_right", 0xc0, -0.0)
161         dist_l_r = 0.285 # Distance between both DWM1000
162         rate = rospy.Rate(10)
163         pos = Position()
164         tf_broadcaster = tf.broadcaster.TransformBroadcaster()
165         rospy.Service('/DWM1000/center', DWM1000Center, handle_center_call)
166
167         while not rospy.is_shutdown() and dwleft.is_alive() and dwright.is_alive():
168                 dist_left = dwleft.distance_valid()
169                 dist_right = dwright.distance_valid()
170                 x = None
171                 y = None
172                 if dist_left == None or dist_right == None:
173                         rospy.logerr_throttle(10, "no valid sensor update %r %r" % (dist_left, dist_right))
174                         # run kalman prediction only
175                         x,y = pos.filter(None, None)
176                 else:
177                         dir = "left" if (dist_left < dist_right) else "right"
178
179                         diff = abs(dist_left - dist_right)
180                         if diff >= dist_l_r:
181                                 # difference to high, correct to maximum
182                                 off = diff - dist_l_r + 0.01
183                                 if dist_left > dist_right:
184                                         dist_left -= off/2
185                                         dist_right += off/2
186                                 else:
187                                         dist_left += off/2
188                                         dist_right -= off/2
189                         rospy.logdebug("%.2f %.2f %.2f %.2f %s", dwleft.distance(), dwright.distance(), dist_left, dist_right, dir)
190
191                         a_r = (-dist_right**2 + dist_left**2 - dist_l_r**2) / (-2*dist_l_r)
192                         x = dist_l_r/2 - a_r
193                         t = dist_right**2 - a_r**2
194                         if t >= 0:
195                                 y = sqrt(t)
196                                 rospy.logdebug("x=%.2f, y=%.2f", x, y)
197                                 # Rotate 90 deg
198                                 x, y = (y, -x)
199
200                                 x, y = pos.filter(x, y)
201
202                                 if VISULAIZE:
203                                         circle_left = plt.Circle((-dist_l_r/2, 0), dwleft.distance, color='red', fill=False)
204                                         circle_right = plt.Circle((dist_l_r/2, 0), dwright.distance, color='green', fill=False)
205                                         plt.gca().add_patch(circle_left)
206                                         plt.gca().add_patch(circle_right)
207                                         plt.grid(True)
208                                         plt.axis('scaled')
209                                         plt.show()
210                         else:
211                                 # No current position, still need up update kalman prediction
212                                 x, y = pos.filter(None, None)
213                 tf_broadcaster.sendTransform((x, y, 0.0), (0, 0, 0, 1), rospy.Time.now(), "uwb_beacon", "base_footprint")
214
215                 rate.sleep()