2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								"""
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								GTSAM Copyright 2010-2020, Georgia Tech Research Corporation,
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								Atlanta, Georgia 30332-0415
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								All Rights Reserved
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								Authors: Frank Dellaert, et al. (see THANKS for the full author list)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								See LICENSE for the license information
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								Track a moving object "Time of Arrival" measurements at 4 microphones.
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								Author: Frank Dellaert
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								"""
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								# pylint: disable=invalid-name, no-name-in-module
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								from gtsam import (LevenbergMarquardtOptimizer, LevenbergMarquardtParams,
							 | 
						
					
						
							
								
									
										
										
										
											2020-08-18 23:02:35 +08:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								                   NonlinearFactorGraph, Point3, Values, noiseModel)
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								from gtsam_unstable import Event, TimeOfArrival, TOAFactor
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								# units
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								MS = 1e-3
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								CM = 1e-2
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								# Instantiate functor with speed of sound value
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								TIME_OF_ARRIVAL = TimeOfArrival(330)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								def define_microphones():
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    """Create microphones."""
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    height = 0.5
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    microphones = []
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    microphones.append(Point3(0, 0, height))
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    microphones.append(Point3(403 * CM, 0, height))
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    microphones.append(Point3(403 * CM, 403 * CM, height))
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    microphones.append(Point3(0, 403 * CM, 2 * height))
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    return microphones
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								def create_trajectory(n):
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    """Create ground truth trajectory."""
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    trajectory = []
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    timeOfEvent = 10
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    # simulate emitting a sound every second while moving on straight line
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    for key in range(n):
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        trajectory.append(
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								            Event(timeOfEvent, 245 * CM + key * 1.0, 201.5 * CM, (212 - 45) * CM))
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        timeOfEvent += 1
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    return trajectory
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								def simulate_one_toa(microphones, event):
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    """Simulate time-of-arrival measurements for a single event."""
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    return [TIME_OF_ARRIVAL.measure(event, microphones[i])
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								            for i in range(len(microphones))]
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								def simulate_toa(microphones, trajectory):
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    """Simulate time-of-arrival measurements for an entire trajectory."""
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    return [simulate_one_toa(microphones, event)
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								            for event in trajectory]
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								def create_graph(microphones, simulatedTOA):
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    """Create factor graph."""
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    graph = NonlinearFactorGraph()
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    # Create a noise model for the TOA error
							 | 
						
					
						
							
								
									
										
										
										
											2020-08-18 23:02:35 +08:00
										 
									 
								 
							 | 
							
								
									
										
									
								
							 | 
							
								
							 | 
							
							
								    model = noiseModel.Isotropic.Sigma(1, 0.5 * MS)
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    K = len(microphones)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    key = 0
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    for toa in simulatedTOA:
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        for i in range(K):
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								            factor = TOAFactor(key, microphones[i], toa[i], model)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								            graph.push_back(factor)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        key += 1
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    return graph
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								def create_initial_estimate(n):
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    """Create initial estimate for n events."""
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    initial = Values()
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    zero = Event()
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    for key in range(n):
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        TOAFactor.InsertEvent(key, zero, initial)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    return initial
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								def toa_example():
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    """Run example with 4 microphones and 5 events in a straight line."""
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    # Create microphones
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    microphones = define_microphones()
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    K = len(microphones)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    for i in range(K):
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        print("mic {} = {}".format(i, microphones[i]))
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    # Create a ground truth trajectory
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    n = 5
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    groundTruth = create_trajectory(n)
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    for event in groundTruth:
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        print(event)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    # Simulate time-of-arrival measurements
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    simulatedTOA = simulate_toa(microphones, groundTruth)
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    for key in range(n):
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								        for i in range(K):
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								            print("z_{}{} = {} ms".format(key, i, simulatedTOA[key][i] / MS))
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    # create factor graph
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    graph = create_graph(microphones, simulatedTOA)
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    print(graph.at(0))
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    # Create initial estimate
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    initial_estimate = create_initial_estimate(n)
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    print(initial_estimate)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    # Optimize using Levenberg-Marquardt optimization.
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    params = LevenbergMarquardtParams()
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    params.setAbsoluteErrorTol(1e-10)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    params.setVerbosityLM("SUMMARY")
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    optimizer = LevenbergMarquardtOptimizer(graph, initial_estimate, params)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    result = optimizer.optimize()
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    print("Final Result:\n", result)
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								
							 | 
						
					
						
							| 
								
							 | 
							
								
							 | 
							
								
							 | 
							
							
								if __name__ == '__main__':
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-25 20:59:58 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    toa_example()
							 | 
						
					
						
							
								
									
										
										
										
											2020-03-19 05:28:40 +08:00
										 
									 
								 
							 | 
							
								
							 | 
							
								
							 | 
							
							
								    print("Example complete")
							 |