so2.py 5.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186
  1. import sympy
  2. import sys
  3. import unittest
  4. import sophus
  5. import functools
  6. class So2:
  7. """ 2 dimensional group of orthogonal matrices with determinant 1 """
  8. def __init__(self, z):
  9. """ internally represented by a unit complex number z """
  10. self.z = z
  11. @staticmethod
  12. def exp(theta):
  13. """ exponential map """
  14. return So2(
  15. sophus.Complex(
  16. sympy.cos(theta),
  17. sympy.sin(theta)))
  18. def log(self):
  19. """ logarithmic map"""
  20. return sympy.atan2(self.z.imag, self.z.real)
  21. def __repr__(self):
  22. return "So2:" + repr(self.z)
  23. @staticmethod
  24. def hat(theta):
  25. return sympy.Matrix([[0, -theta],
  26. [theta, 0]])
  27. def matrix(self):
  28. """ returns matrix representation """
  29. return sympy.Matrix([
  30. [self.z.real, -self.z.imag],
  31. [self.z.imag, self.z.real]])
  32. def __mul__(self, right):
  33. """ left-multiplication
  34. either rotation concatenation or point-transform """
  35. if isinstance(right, sympy.Matrix):
  36. assert right.shape == (2, 1), right.shape
  37. return self.matrix() * right
  38. elif isinstance(right, So2):
  39. return So2(self.z * right.z)
  40. assert False, "unsupported type: {0}".format(type(right))
  41. def __getitem__(self, key):
  42. return self.z[key]
  43. @staticmethod
  44. def calc_Dx_exp_x(x):
  45. return sympy.Matrix(2, 1, lambda r, c:
  46. sympy.diff(So2.exp(x)[r], x))
  47. @staticmethod
  48. def Dx_exp_x_at_0():
  49. return sympy.Matrix([0, 1])
  50. @staticmethod
  51. def calc_Dx_exp_x_at_0(x):
  52. return So2.calc_Dx_exp_x(x).limit(x, 0)
  53. def calc_Dx_this_mul_exp_x_at_0(self, x):
  54. return sympy.Matrix(2, 1, lambda r, c:
  55. sympy.diff((self * So2.exp(x))[r], x))\
  56. .limit(x, 0)
  57. @staticmethod
  58. def Dxi_x_matrix(x, i):
  59. if i == 0:
  60. return sympy.Matrix([[1, 0],
  61. [0, 1]])
  62. if i == 1:
  63. return sympy.Matrix([[0, -1],
  64. [1, 0]])
  65. @staticmethod
  66. def calc_Dxi_x_matrix(x, i):
  67. return sympy.Matrix(2, 2, lambda r, c:
  68. sympy.diff(x.matrix()[r, c], x[i]))
  69. @staticmethod
  70. def Dx_exp_x_matrix(x):
  71. R = So2.exp(x)
  72. Dx_exp_x = So2.calc_Dx_exp_x(x)
  73. l = [Dx_exp_x[j] * So2.Dxi_x_matrix(R, j) for j in [0, 1]]
  74. return functools.reduce((lambda a, b: a + b), l)
  75. @staticmethod
  76. def calc_Dx_exp_x_matrix(x):
  77. return sympy.Matrix(2, 2, lambda r, c:
  78. sympy.diff(So2.exp(x).matrix()[r, c], x))
  79. @staticmethod
  80. def Dx_exp_x_matrix_at_0():
  81. return So2.hat(1)
  82. @staticmethod
  83. def calc_Dx_exp_x_matrix_at_0(x):
  84. return sympy.Matrix(2, 2, lambda r, c:
  85. sympy.diff(So2.exp(x).matrix()[r, c], x)
  86. ).limit(x, 0)
  87. class TestSo2(unittest.TestCase):
  88. def setUp(self):
  89. self.theta = sympy.symbols(
  90. 'theta', real=True)
  91. x, y = sympy.symbols('c[0] c[1]', real=True)
  92. p0, p1 = sympy.symbols('p0 p1', real=True)
  93. self.a = So2(sophus.Complex(x, y))
  94. self.p = sophus.Vector2(p0, p1)
  95. def test_exp_log(self):
  96. for theta in [0., 0.5, 0.1]:
  97. w = So2.exp(theta).log()
  98. self.assertAlmostEqual(theta, w)
  99. def test_matrix(self):
  100. R_foo_bar = So2.exp(self.theta)
  101. Rmat_foo_bar = R_foo_bar.matrix()
  102. point_bar = self.p
  103. p1_foo = R_foo_bar * point_bar
  104. p2_foo = Rmat_foo_bar * point_bar
  105. self.assertEqual(sympy.simplify(p1_foo - p2_foo),
  106. sophus.ZeroVector2())
  107. def test_derivatives(self):
  108. self.assertEqual(sympy.simplify(So2.calc_Dx_exp_x_at_0(self.theta) -
  109. So2.Dx_exp_x_at_0()),
  110. sympy.Matrix.zeros(2, 1))
  111. for i in [0, 1]:
  112. self.assertEqual(sympy.simplify(So2.calc_Dxi_x_matrix(self.a, i) -
  113. So2.Dxi_x_matrix(self.a, i)),
  114. sympy.Matrix.zeros(2, 2))
  115. self.assertEqual(sympy.simplify(
  116. So2.Dx_exp_x_matrix(self.theta) -
  117. So2.calc_Dx_exp_x_matrix(self.theta)),
  118. sympy.Matrix.zeros(2, 2))
  119. self.assertEqual(sympy.simplify(
  120. So2.Dx_exp_x_matrix_at_0() -
  121. So2.calc_Dx_exp_x_matrix_at_0(self.theta)),
  122. sympy.Matrix.zeros(2, 2))
  123. def test_codegen(self):
  124. stream = sophus.cse_codegen(So2.calc_Dx_exp_x(self.theta))
  125. filename = "cpp_gencode/So2_Dx_exp_x.cpp"
  126. # set to true to generate codegen files
  127. if False:
  128. file = open(filename, "w")
  129. for line in stream:
  130. file.write(line)
  131. file.close()
  132. else:
  133. file = open(filename, "r")
  134. file_lines = file.readlines()
  135. for i, line in enumerate(stream):
  136. self.assertEqual(line, file_lines[i])
  137. file.close()
  138. stream.close
  139. stream = sophus.cse_codegen(
  140. self.a.calc_Dx_this_mul_exp_x_at_0(self.theta))
  141. filename = "cpp_gencode/So2_Dx_this_mul_exp_x_at_0.cpp"
  142. # set to true to generate codegen files
  143. if False:
  144. file = open(filename, "w")
  145. for line in stream:
  146. file.write(line)
  147. file.close()
  148. else:
  149. file = open(filename, "r")
  150. file_lines = file.readlines()
  151. for i, line in enumerate(stream):
  152. self.assertEqual(line, file_lines[i])
  153. file.close()
  154. stream.close
  155. if __name__ == '__main__':
  156. unittest.main()