Kalman.cpp 37 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873
  1. #include "Kalman.h"
  2. #include <iostream>
  3. #include <vector>
  4. //---------------------------------------------------------------------------
  5. TKalmanFilter::TKalmanFilter(
  6. tracking::KalmanType type,
  7. bool useAcceleration,
  8. track_t deltaTime, // time increment (lower values makes target more "massive")
  9. track_t accelNoiseMag
  10. )
  11. :
  12. m_type(type),
  13. m_initialized(false),
  14. m_deltaTime(deltaTime),
  15. m_deltaTimeMin(deltaTime),
  16. m_deltaTimeMax(2 * deltaTime),
  17. m_lastDist(0),
  18. m_accelNoiseMag(accelNoiseMag),
  19. m_useAcceleration(useAcceleration)
  20. {
  21. m_deltaStep = (m_deltaTimeMax - m_deltaTimeMin) / m_deltaStepsCount;
  22. }
  23. //---------------------------------------------------------------------------
  24. void TKalmanFilter::CreateLinear(cv::Point3f xy0, cv::Point3f xyv0)
  25. {
  26. // We don't know acceleration, so, assume it to process noise.
  27. // But we can guess, the range of acceleration values thich can be achieved by tracked object.
  28. // Process noise. (standard deviation of acceleration: m/s^2)
  29. // shows, woh much target can accelerate.
  30. // 4 state variables, 2 measurements
  31. m_linearKalman.init(4, 2, 0, El_t);
  32. // Transition cv::Matrix
  33. m_linearKalman.transitionMatrix = (cv::Mat_<track_t>(4, 4) <<
  34. 1, 0, m_deltaTime, 0,
  35. 0, 1, 0, m_deltaTime,
  36. 0, 0, 1, 0,
  37. 0, 0, 0, 1);
  38. // init...
  39. m_lastPointResult = xy0;
  40. m_linearKalman.statePre.at<track_t>(0) = xy0.x; // x
  41. m_linearKalman.statePre.at<track_t>(1) = xy0.y; // y
  42. m_linearKalman.statePre.at<track_t>(2) = xyv0.x; // vx
  43. m_linearKalman.statePre.at<track_t>(3) = xyv0.y; // vy
  44. m_linearKalman.statePost.at<track_t>(0) = xy0.x;
  45. m_linearKalman.statePost.at<track_t>(1) = xy0.y;
  46. m_linearKalman.statePost.at<track_t>(2) = xyv0.x;
  47. m_linearKalman.statePost.at<track_t>(3) = xyv0.y;
  48. cv::setIdentity(m_linearKalman.measurementMatrix);
  49. m_linearKalman.processNoiseCov = (cv::Mat_<track_t>(4, 4) <<
  50. pow(m_deltaTime,4.0)/4.0 ,0 ,pow(m_deltaTime,3.0)/2.0 ,0,
  51. 0 ,pow(m_deltaTime,4.0)/4.0 ,0 ,pow(m_deltaTime,3.0)/2.0,
  52. pow(m_deltaTime,3.0)/2.0 ,0 ,pow(m_deltaTime,2.0) ,0,
  53. 0 ,pow(m_deltaTime,3.0)/2.0 ,0 ,pow(m_deltaTime,2.0));
  54. m_linearKalman.processNoiseCov *= m_accelNoiseMag;
  55. cv::setIdentity(m_linearKalman.measurementNoiseCov, cv::Scalar::all(0.1));
  56. cv::setIdentity(m_linearKalman.errorCovPost, cv::Scalar::all(.1));
  57. m_initialized = true;
  58. }
  59. //---------------------------------------------------------------------------
  60. void TKalmanFilter::CreateLinear(cv::Rect_<track_t> rect0, Point_t rectv0)
  61. {
  62. // We don't know acceleration, so, assume it to process noise.
  63. // But we can guess, the range of acceleration values thich can be achieved by tracked object.
  64. // Process noise. (standard deviation of acceleration: m/s^2)
  65. // shows, woh much target can accelerate.
  66. // 8 state variables (x, y, vx, vy, width, height, vw, vh), 4 measurements (x, y, width, height)
  67. m_linearKalman.init(8, 4, 0, El_t);
  68. // Transition cv::Matrix
  69. m_linearKalman.transitionMatrix = (cv::Mat_<track_t>(8, 8) <<
  70. 1, 0, 0, 0, m_deltaTime, 0, 0, 0,
  71. 0, 1, 0, 0, 0, m_deltaTime, 0, 0,
  72. 0, 0, 1, 0, 0, 0, m_deltaTime, 0,
  73. 0, 0, 0, 1, 0, 0, 0, m_deltaTime,
  74. 0, 0, 0, 0, 1, 0, 0, 0,
  75. 0, 0, 0, 0, 0, 1, 0, 0,
  76. 0, 0, 0, 0, 0, 0, 1, 0,
  77. 0, 0, 0, 0, 0, 0, 0, 1);
  78. // init...
  79. m_linearKalman.statePre.at<track_t>(0) = rect0.x; // x
  80. m_linearKalman.statePre.at<track_t>(1) = rect0.y; // y
  81. m_linearKalman.statePre.at<track_t>(2) = rect0.width; // width
  82. m_linearKalman.statePre.at<track_t>(3) = rect0.height; // height
  83. m_linearKalman.statePre.at<track_t>(4) = rectv0.x; // dx
  84. m_linearKalman.statePre.at<track_t>(5) = rectv0.y; // dy
  85. m_linearKalman.statePre.at<track_t>(6) = 0; // dw
  86. m_linearKalman.statePre.at<track_t>(7) = 0; // dh
  87. m_linearKalman.statePost.at<track_t>(0) = rect0.x;
  88. m_linearKalman.statePost.at<track_t>(1) = rect0.y;
  89. m_linearKalman.statePost.at<track_t>(2) = rect0.width;
  90. m_linearKalman.statePost.at<track_t>(3) = rect0.height;
  91. m_linearKalman.statePost.at<track_t>(4) = rectv0.x;
  92. m_linearKalman.statePost.at<track_t>(5) = rectv0.y;
  93. m_linearKalman.statePost.at<track_t>(6) = 0;
  94. m_linearKalman.statePost.at<track_t>(7) = 0;
  95. cv::setIdentity(m_linearKalman.measurementMatrix);
  96. track_t n1 = pow(m_deltaTime, 4.f) / 4.f;
  97. track_t n2 = pow(m_deltaTime, 3.f) / 2.f;
  98. track_t n3 = pow(m_deltaTime, 2.f);
  99. m_linearKalman.processNoiseCov = (cv::Mat_<track_t>(8, 8) <<
  100. n1, 0, 0, 0, n2, 0, 0, 0,
  101. 0, n1, 0, 0, 0, n2, 0, 0,
  102. 0, 0, n1, 0, 0, 0, n2, 0,
  103. 0, 0, 0, n1, 0, 0, 0, n2,
  104. n2, 0, 0, 0, n3, 0, 0, 0,
  105. 0, n2, 0, 0, 0, n3, 0, 0,
  106. 0, 0, n2, 0, 0, 0, n3, 0,
  107. 0, 0, 0, n2, 0, 0, 0, n3);
  108. m_linearKalman.processNoiseCov *= m_accelNoiseMag;
  109. cv::setIdentity(m_linearKalman.measurementNoiseCov, cv::Scalar::all(0.1));
  110. cv::setIdentity(m_linearKalman.errorCovPost, cv::Scalar::all(.1));
  111. m_initialized = true;
  112. }
  113. //---------------------------------------------------------------------------
  114. void TKalmanFilter::CreateLinear3D(Rect3D rect0, cv::Point3f rectv0)
  115. {
  116. // We don't know acceleration, so, assume it to process noise.
  117. // But we can guess, the range of acceleration values thich can be achieved by tracked object.
  118. // Process noise. (standard deviation of acceleration: m/s^2)
  119. // shows, woh much target can accelerate.
  120. // 14 state variables (x, y, z, width, height, length, yaw, d...), 7 measurements (x, y, z, width, height, length, yaw)
  121. m_linearKalman.init(14, 7, 0, El_t);
  122. // Transition cv::Matrix
  123. m_linearKalman.transitionMatrix = (cv::Mat_<track_t>(14, 14) <<
  124. 1, 0, 0, 0, 0, 0, 0, m_deltaTime, 0, 0, 0, 0, 0, 0,
  125. 0, 1, 0, 0, 0, 0, 0, 0, m_deltaTime, 0, 0, 0, 0, 0,
  126. 0, 0, 1, 0, 0, 0, 0, 0, 0, m_deltaTime, 0, 0, 0, 0,
  127. 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, m_deltaTime, 0, 0, 0,
  128. 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, m_deltaTime, 0, 0,
  129. 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, m_deltaTime, 0,
  130. 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, m_deltaTime,
  131. 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
  132. 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
  133. 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
  134. 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
  135. 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
  136. 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
  137. 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1);
  138. // init...
  139. m_linearKalman.statePre.at<track_t>(0) = rect0.center.x; // x
  140. m_linearKalman.statePre.at<track_t>(1) = rect0.center.y; // y
  141. m_linearKalman.statePre.at<track_t>(2) = rect0.center.z; // z
  142. m_linearKalman.statePre.at<track_t>(3) = rect0.size.width; // width
  143. m_linearKalman.statePre.at<track_t>(4) = rect0.size.height; // height
  144. m_linearKalman.statePre.at<track_t>(5) = rect0.size.length; // length
  145. m_linearKalman.statePre.at<track_t>(6) = rect0.yaw; // yaw
  146. m_linearKalman.statePre.at<track_t>(7) = rectv0.x; // dx
  147. m_linearKalman.statePre.at<track_t>(8) = rectv0.y; // dy
  148. m_linearKalman.statePre.at<track_t>(9) = 0; // dz
  149. m_linearKalman.statePre.at<track_t>(10) = 0; // dw
  150. m_linearKalman.statePre.at<track_t>(11) = 0; // dh
  151. m_linearKalman.statePre.at<track_t>(12) = 0; // dl
  152. m_linearKalman.statePre.at<track_t>(13) = 0; // dyaw
  153. m_linearKalman.statePost.at<track_t>(0) = rect0.center.x;
  154. m_linearKalman.statePost.at<track_t>(1) = rect0.center.y;
  155. m_linearKalman.statePost.at<track_t>(2) = rect0.center.z;
  156. m_linearKalman.statePost.at<track_t>(3) = rect0.size.width;
  157. m_linearKalman.statePost.at<track_t>(4) = rect0.size.height;
  158. m_linearKalman.statePost.at<track_t>(5) = rect0.size.length;
  159. m_linearKalman.statePost.at<track_t>(6) = rect0.yaw;
  160. m_linearKalman.statePost.at<track_t>(7) = rectv0.x;
  161. m_linearKalman.statePost.at<track_t>(8) = rectv0.y;
  162. m_linearKalman.statePost.at<track_t>(9) = 0;
  163. m_linearKalman.statePost.at<track_t>(10) = 0;
  164. m_linearKalman.statePost.at<track_t>(11) = 0;
  165. m_linearKalman.statePost.at<track_t>(12) = 0;
  166. m_linearKalman.statePost.at<track_t>(13) = 0;
  167. cv::setIdentity(m_linearKalman.measurementMatrix);
  168. track_t n1 = pow(m_deltaTime, 4.f) / 4.f;
  169. track_t n2 = pow(m_deltaTime, 3.f) / 2.f;
  170. track_t n3 = pow(m_deltaTime, 2.f);
  171. m_linearKalman.processNoiseCov = (cv::Mat_<track_t>(14, 14) <<
  172. n1, 0, 0, 0, 0, 0, 0, n2, 0, 0, 0, 0, 0, 0,
  173. 0, n1, 0, 0, 0, 0, 0, 0, n2, 0, 0, 0, 0, 0,
  174. 0, 0, n1, 0, 0, 0, 0, 0, 0, n2, 0, 0, 0, 0,
  175. 0, 0, 0, n1, 0, 0, 0, 0, 0, 0, n2, 0, 0, 0,
  176. 0, 0, 0, 0, n1, 0, 0, 0, 0, 0, 0, n2, 0, 0,
  177. 0, 0, 0, 0, 0, n1, 0, 0, 0, 0, 0, 0, n2, 0,
  178. 0, 0, 0, 0, 0, 0, n1, 0, 0, 0, 0, 0, 0, n2,
  179. n2, 0, 0, 0, 0, 0, 0, n3, 0, 0, 0, 0, 0, 0,
  180. 0, n2, 0, 0, 0, 0, 0, 0, n3, 0, 0, 0, 0, 0,
  181. 0, 0, n2, 0, 0, 0, 0, 0, 0, n3, 0, 0, 0, 0,
  182. 0, 0, 0, n2, 0, 0, 0, 0, 0, 0, n3, 0, 0, 0,
  183. 0, 0, 0, 0, n2, 0, 0, 0, 0, 0, 0, n3, 0, 0,
  184. 0, 0, 0, 0, 0, n2, 0, 0, 0, 0, 0, 0, n3, 0,
  185. 0, 0, 0, 0, 0, 0, n2, 0, 0, 0, 0, 0, 0, n3);
  186. m_linearKalman.processNoiseCov *= m_accelNoiseMag;
  187. cv::setIdentity(m_linearKalman.measurementNoiseCov, cv::Scalar::all(0.1));
  188. cv::setIdentity(m_linearKalman.errorCovPost, cv::Scalar::all(.1));
  189. m_initialized = true;
  190. }
  191. //---------------------------------------------------------------------------
  192. void TKalmanFilter::CreateLinearAcceleration(cv::Point3f xy0, cv::Point3f xyv0)
  193. {
  194. // 6 state variables, 2 measurements
  195. m_linearKalman.init(6, 2, 0, El_t);
  196. // Transition cv::Matrix
  197. const track_t dt = m_deltaTime;
  198. const track_t dt2 = 0.5f * m_deltaTime * m_deltaTime;
  199. m_linearKalman.transitionMatrix = (cv::Mat_<track_t>(6, 6) <<
  200. 1, 0, dt, 0, dt2, 0,
  201. 0, 1, 0, dt, 0, dt2,
  202. 0, 0, 1, 0, dt, 0,
  203. 0, 0, 0, 1, 0, dt,
  204. 0, 0, 0, 0, 1, 0,
  205. 0, 0, 0, 0, 0, 1);
  206. // init...
  207. m_lastPointResult = xy0;
  208. m_linearKalman.statePre.at<track_t>(0) = xy0.x; // x
  209. m_linearKalman.statePre.at<track_t>(1) = xy0.y; // y
  210. m_linearKalman.statePre.at<track_t>(2) = xyv0.x; // vx
  211. m_linearKalman.statePre.at<track_t>(3) = xyv0.y; // vy
  212. m_linearKalman.statePre.at<track_t>(4) = 0; // ax
  213. m_linearKalman.statePre.at<track_t>(5) = 0; // ay
  214. m_linearKalman.statePost.at<track_t>(0) = xy0.x;
  215. m_linearKalman.statePost.at<track_t>(1) = xy0.y;
  216. m_linearKalman.statePost.at<track_t>(2) = xyv0.x;
  217. m_linearKalman.statePost.at<track_t>(3) = xyv0.y;
  218. m_linearKalman.statePost.at<track_t>(4) = 0;
  219. m_linearKalman.statePost.at<track_t>(5) = 0;
  220. cv::setIdentity(m_linearKalman.measurementMatrix);
  221. track_t n1 = pow(m_deltaTime, 4.f) / 4.f;
  222. track_t n2 = pow(m_deltaTime, 3.f) / 2.f;
  223. track_t n3 = pow(m_deltaTime, 2.f);
  224. m_linearKalman.processNoiseCov = (cv::Mat_<track_t>(6, 6) <<
  225. n1, 0, n2, 0, n2, 0,
  226. 0, n1, 0, n2, 0, n2,
  227. n2, 0, n3, 0, n3, 0,
  228. 0, n2, 0, n3, 0, n3,
  229. 0, 0, n2, 0, n3, 0,
  230. 0, 0, 0, n2, 0, n3);
  231. m_linearKalman.processNoiseCov *= m_accelNoiseMag;
  232. cv::setIdentity(m_linearKalman.measurementNoiseCov, cv::Scalar::all(0.1));
  233. cv::setIdentity(m_linearKalman.errorCovPost, cv::Scalar::all(.1));
  234. m_initialized = true;
  235. }
  236. //---------------------------------------------------------------------------
  237. void TKalmanFilter::CreateLinearAcceleration(cv::Rect_<track_t> rect0, Point_t rectv0)
  238. {
  239. // 12 state variables (x, y, vx, vy, ax, ay, width, height, vw, vh, aw, ah), 4 measurements (x, y, width, height)
  240. m_linearKalman.init(12, 4, 0, El_t);
  241. // Transition cv::Matrix
  242. const track_t dt = m_deltaTime;
  243. const track_t dt2 = 0.5f * m_deltaTime * m_deltaTime;
  244. m_linearKalman.transitionMatrix = (cv::Mat_<track_t>(12, 12) <<
  245. 1, 0, 0, 0, dt, 0, 0, 0, dt2, 0, dt2, 0,
  246. 0, 1, 0, 0, 0, dt, 0, 0, 0, dt2, 0, dt2,
  247. 0, 0, 1, 0, 0, 0, dt, 0, 0, 0, dt2, 0,
  248. 0, 0, 0, 1, 0, 0, 0, dt, 0, 0, 0, dt2,
  249. 0, 0, 0, 0, 1, 0, 0, 0, dt, 0, 0, 0,
  250. 0, 0, 0, 0, 0, 1, 0, 0, 0, dt, 0, 0,
  251. 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, dt, 0,
  252. 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, dt,
  253. 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
  254. 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
  255. 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1);
  256. // init...
  257. m_linearKalman.statePre.at<track_t>(0) = rect0.x; // x
  258. m_linearKalman.statePre.at<track_t>(1) = rect0.y; // y
  259. m_linearKalman.statePre.at<track_t>(2) = rect0.width; // width
  260. m_linearKalman.statePre.at<track_t>(3) = rect0.height; // height
  261. m_linearKalman.statePre.at<track_t>(4) = rectv0.x; // dx
  262. m_linearKalman.statePre.at<track_t>(5) = rectv0.y; // dy
  263. m_linearKalman.statePre.at<track_t>(6) = 0; // dw
  264. m_linearKalman.statePre.at<track_t>(7) = 0; // dh
  265. m_linearKalman.statePre.at<track_t>(8) = 0; // ax
  266. m_linearKalman.statePre.at<track_t>(9) = 0; // ay
  267. m_linearKalman.statePre.at<track_t>(10) = 0; // aw
  268. m_linearKalman.statePre.at<track_t>(11) = 0; // ah
  269. m_linearKalman.statePost.at<track_t>(0) = rect0.x;
  270. m_linearKalman.statePost.at<track_t>(1) = rect0.y;
  271. m_linearKalman.statePost.at<track_t>(2) = rect0.width;
  272. m_linearKalman.statePost.at<track_t>(3) = rect0.height;
  273. m_linearKalman.statePost.at<track_t>(4) = rectv0.x;
  274. m_linearKalman.statePost.at<track_t>(5) = rectv0.y;
  275. m_linearKalman.statePost.at<track_t>(6) = 0;
  276. m_linearKalman.statePost.at<track_t>(7) = 0;
  277. m_linearKalman.statePost.at<track_t>(8) = 0;
  278. m_linearKalman.statePost.at<track_t>(9) = 0;
  279. m_linearKalman.statePost.at<track_t>(10) = 0;
  280. m_linearKalman.statePost.at<track_t>(11) = 0;
  281. cv::setIdentity(m_linearKalman.measurementMatrix);
  282. track_t n1 = pow(m_deltaTime, 4.f) / 4.f;
  283. track_t n2 = pow(m_deltaTime, 3.f) / 2.f;
  284. track_t n3 = pow(m_deltaTime, 2.f);
  285. m_linearKalman.processNoiseCov = (cv::Mat_<track_t>(12, 12) <<
  286. n1, 0, 0, 0, n2, 0, 0, 0, n2, 0, n2,
  287. 0, n1, 0, 0, 0, n2, 0, 0, 0, n2, 0, n2,
  288. 0, 0, n1, 0, 0, 0, n2, 0, 0, 0, n2, 0,
  289. 0, 0, 0, n1, 0, 0, 0, n2, 0, 0, 0, n2,
  290. n2, 0, 0, 0, n3, 0, 0, 0, n3, 0, n3, 0,
  291. 0, n2, 0, 0, 0, n3, 0, 0, 0, n3, 0, n3,
  292. 0, 0, n2, 0, 0, 0, n3, 0, 0, 0, n3, 0,
  293. 0, 0, 0, n2, 0, 0, 0, n3, 0, 0, 0, n3,
  294. n2, 0, 0, 0, n3, 0, 0, 0, n3, 0, 0, 0,
  295. 0, n2, 0, 0, 0, n3, 0, 0, 0, n3, 0, 0,
  296. 0, 0, n2, 0, 0, 0, n3, 0, 0, 0, n3, 0,
  297. 0, 0, 0, n2, 0, 0, 0, n3, 0, 0, 0, n3);
  298. m_linearKalman.processNoiseCov *= m_accelNoiseMag;
  299. cv::setIdentity(m_linearKalman.measurementNoiseCov, cv::Scalar::all(0.1));
  300. cv::setIdentity(m_linearKalman.errorCovPost, cv::Scalar::all(.1));
  301. m_initialized = true;
  302. }
  303. //---------------------------------------------------------------------------
  304. void TKalmanFilter::CreateLinearAcceleration3D(Rect3D rect0, cv::Point3f rectv0)
  305. {
  306. // 12 state variables (x, y, vx, vy, ax, ay, width, height, vw, vh, aw, ah), 4 measurements (x, y, width, height)
  307. m_linearKalman.init(12, 4, 0, El_t);
  308. // Transition cv::Matrix
  309. const track_t dt = m_deltaTime;
  310. const track_t dt2 = 0.5f * m_deltaTime * m_deltaTime;
  311. m_linearKalman.transitionMatrix = (cv::Mat_<track_t>(12, 12) <<
  312. 1, 0, 0, 0, dt, 0, 0, 0, dt2, 0, dt2, 0,
  313. 0, 1, 0, 0, 0, dt, 0, 0, 0, dt2, 0, dt2,
  314. 0, 0, 1, 0, 0, 0, dt, 0, 0, 0, dt2, 0,
  315. 0, 0, 0, 1, 0, 0, 0, dt, 0, 0, 0, dt2,
  316. 0, 0, 0, 0, 1, 0, 0, 0, dt, 0, 0, 0,
  317. 0, 0, 0, 0, 0, 1, 0, 0, 0, dt, 0, 0,
  318. 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, dt, 0,
  319. 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, dt,
  320. 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
  321. 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
  322. 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1);
  323. // init...
  324. m_linearKalman.statePre.at<track_t>(0) = rect0.center.x; // x
  325. m_linearKalman.statePre.at<track_t>(1) = rect0.center.y; // y
  326. m_linearKalman.statePre.at<track_t>(2) = rect0.size.width; // width
  327. m_linearKalman.statePre.at<track_t>(3) = rect0.size.height; // height
  328. m_linearKalman.statePre.at<track_t>(4) = rectv0.x; // dx
  329. m_linearKalman.statePre.at<track_t>(5) = rectv0.y; // dy
  330. m_linearKalman.statePre.at<track_t>(6) = 0; // dw
  331. m_linearKalman.statePre.at<track_t>(7) = 0; // dh
  332. m_linearKalman.statePre.at<track_t>(8) = 0; // ax
  333. m_linearKalman.statePre.at<track_t>(9) = 0; // ay
  334. m_linearKalman.statePre.at<track_t>(10) = 0; // aw
  335. m_linearKalman.statePre.at<track_t>(11) = 0; // ah
  336. m_linearKalman.statePost.at<track_t>(0) = rect0.center.x;
  337. m_linearKalman.statePost.at<track_t>(1) = rect0.center.y;
  338. m_linearKalman.statePost.at<track_t>(2) = rect0.size.width;
  339. m_linearKalman.statePost.at<track_t>(3) = rect0.size.height;
  340. m_linearKalman.statePost.at<track_t>(4) = rectv0.x;
  341. m_linearKalman.statePost.at<track_t>(5) = rectv0.y;
  342. m_linearKalman.statePost.at<track_t>(6) = 0;
  343. m_linearKalman.statePost.at<track_t>(7) = 0;
  344. m_linearKalman.statePost.at<track_t>(8) = 0;
  345. m_linearKalman.statePost.at<track_t>(9) = 0;
  346. m_linearKalman.statePost.at<track_t>(10) = 0;
  347. m_linearKalman.statePost.at<track_t>(11) = 0;
  348. cv::setIdentity(m_linearKalman.measurementMatrix);
  349. track_t n1 = pow(m_deltaTime, 4.f) / 4.f;
  350. track_t n2 = pow(m_deltaTime, 3.f) / 2.f;
  351. track_t n3 = pow(m_deltaTime, 2.f);
  352. m_linearKalman.processNoiseCov = (cv::Mat_<track_t>(12, 12) <<
  353. n1, 0, 0, 0, n2, 0, 0, 0, n2, 0, n2,
  354. 0, n1, 0, 0, 0, n2, 0, 0, 0, n2, 0, n2,
  355. 0, 0, n1, 0, 0, 0, n2, 0, 0, 0, n2, 0,
  356. 0, 0, 0, n1, 0, 0, 0, n2, 0, 0, 0, n2,
  357. n2, 0, 0, 0, n3, 0, 0, 0, n3, 0, n3, 0,
  358. 0, n2, 0, 0, 0, n3, 0, 0, 0, n3, 0, n3,
  359. 0, 0, n2, 0, 0, 0, n3, 0, 0, 0, n3, 0,
  360. 0, 0, 0, n2, 0, 0, 0, n3, 0, 0, 0, n3,
  361. n2, 0, 0, 0, n3, 0, 0, 0, n3, 0, 0, 0,
  362. 0, n2, 0, 0, 0, n3, 0, 0, 0, n3, 0, 0,
  363. 0, 0, n2, 0, 0, 0, n3, 0, 0, 0, n3, 0,
  364. 0, 0, 0, n2, 0, 0, 0, n3, 0, 0, 0, n3);
  365. m_linearKalman.processNoiseCov *= m_accelNoiseMag;
  366. cv::setIdentity(m_linearKalman.measurementNoiseCov, cv::Scalar::all(0.1));
  367. cv::setIdentity(m_linearKalman.errorCovPost, cv::Scalar::all(.1));
  368. m_initialized = true;
  369. }
  370. //---------------------------------------------------------------------------
  371. cv::Point3f TKalmanFilter::GetPointPrediction()
  372. {
  373. if (m_initialized)
  374. {
  375. cv::Mat prediction;
  376. switch (m_type)
  377. {
  378. case tracking::KalmanLinear:
  379. prediction = m_linearKalman.predict();
  380. break;
  381. }
  382. m_lastPointResult = cv::Point3f(prediction.at<track_t>(0), prediction.at<track_t>(1), prediction.at<track_t>(2));
  383. }
  384. else
  385. {
  386. }
  387. return m_lastPointResult;
  388. }
  389. //---------------------------------------------------------------------------
  390. cv::Point3f TKalmanFilter::Update(cv::Point3f pt, bool dataCorrect)
  391. {
  392. if (!m_initialized)
  393. {
  394. if (m_initialPoints.size() < MIN_INIT_VALS)
  395. {
  396. if (dataCorrect)
  397. {
  398. m_initialPoints.push_back(pt);
  399. m_lastPointResult = pt;
  400. }
  401. }
  402. if (m_initialPoints.size() == MIN_INIT_VALS)
  403. {
  404. track_t kx = 0;
  405. track_t bx = 0;
  406. track_t ky = 0;
  407. track_t by = 0;
  408. get_lin_regress_params(m_initialPoints, 0, MIN_INIT_VALS, kx, bx, ky, by);//predict p,v
  409. cv::Point3f xy0(kx * (MIN_INIT_VALS - 1) + bx, ky * (MIN_INIT_VALS - 1) + by, m_lastPointResult.z);
  410. cv::Point3f xyv0(kx, ky,0);
  411. switch (m_type)
  412. {
  413. case tracking::KalmanLinear:
  414. if (m_useAcceleration)
  415. CreateLinearAcceleration(xy0, xyv0);
  416. else
  417. CreateLinear(xy0, xyv0);
  418. break;
  419. }
  420. m_lastDist = 0;
  421. }
  422. }
  423. if (m_initialized)
  424. {
  425. cv::Mat measurement(2, 1, Mat_t(1));
  426. if (!dataCorrect)
  427. {
  428. measurement.at<track_t>(0) = m_lastPointResult.x; //update using prediction
  429. measurement.at<track_t>(1) = m_lastPointResult.y;
  430. }
  431. else
  432. {
  433. measurement.at<track_t>(0) = pt.x; //update using measurements
  434. measurement.at<track_t>(1) = pt.y;
  435. }
  436. // Correction
  437. cv::Mat estimated;
  438. switch (m_type)
  439. {
  440. case tracking::KalmanLinear:
  441. {
  442. estimated = m_linearKalman.correct(measurement);
  443. // Inertia correction
  444. if (!m_useAcceleration)
  445. {
  446. track_t currDist = sqrtf(sqr(estimated.at<track_t>(0) - pt.x) + sqr(estimated.at<track_t>(1) - pt.y));
  447. if (currDist > m_lastDist)
  448. {
  449. m_deltaTime = std::min(m_deltaTime + m_deltaStep, m_deltaTimeMax);
  450. }
  451. else
  452. {
  453. m_deltaTime = std::max(m_deltaTime - m_deltaStep, m_deltaTimeMin);
  454. }
  455. m_lastDist = currDist;
  456. m_linearKalman.transitionMatrix.at<track_t>(0, 2) = m_deltaTime;
  457. m_linearKalman.transitionMatrix.at<track_t>(1, 3) = m_deltaTime;
  458. }
  459. break;
  460. }
  461. }
  462. m_lastPointResult.x = estimated.at<track_t>(0); //update using measurements
  463. m_lastPointResult.y = estimated.at<track_t>(1);
  464. }
  465. else
  466. {
  467. if (dataCorrect)
  468. {
  469. m_lastPointResult = pt;
  470. }
  471. }
  472. return m_lastPointResult;
  473. }
  474. //---------------------------------------------------------------------------
  475. cv::Rect TKalmanFilter::GetRectPrediction()
  476. {
  477. if (m_initialized)
  478. {
  479. cv::Mat prediction;
  480. switch (m_type)
  481. {
  482. case tracking::KalmanLinear:
  483. prediction = m_linearKalman.predict();
  484. break;
  485. }
  486. m_lastRectResult = cv::Rect_<track_t>(prediction.at<track_t>(0), prediction.at<track_t>(1), prediction.at<track_t>(2), prediction.at<track_t>(3));
  487. }
  488. else
  489. {
  490. }
  491. return cv::Rect(static_cast<int>(m_lastRectResult.x), static_cast<int>(m_lastRectResult.y), static_cast<int>(m_lastRectResult.width), static_cast<int>(m_lastRectResult.height));
  492. }
  493. //---------------------------------------------------------------------------
  494. cv::Rect TKalmanFilter::Update(cv::Rect rect, bool dataCorrect)
  495. {
  496. if (!m_initialized)
  497. {
  498. if (m_initialRects.size() < MIN_INIT_VALS)
  499. {
  500. if (dataCorrect)
  501. {
  502. m_initialRects.push_back(rect);
  503. m_lastRectResult.x = static_cast<track_t>(rect.x);
  504. m_lastRectResult.y = static_cast<track_t>(rect.y);
  505. m_lastRectResult.width = static_cast<track_t>(rect.width);
  506. m_lastRectResult.height = static_cast<track_t>(rect.height);
  507. }
  508. }
  509. if (m_initialRects.size() == MIN_INIT_VALS)
  510. {
  511. std::vector<Point_t> initialPoints;
  512. Point_t averageSize(0, 0);
  513. for (const auto& r : m_initialRects)
  514. {
  515. initialPoints.emplace_back(static_cast<track_t>(r.x), static_cast<track_t>(r.y));
  516. averageSize.x += r.width;
  517. averageSize.y += r.height;
  518. }
  519. averageSize.x /= MIN_INIT_VALS;
  520. averageSize.y /= MIN_INIT_VALS;
  521. track_t kx = 0;
  522. track_t bx = 0;
  523. track_t ky = 0;
  524. track_t by = 0;
  525. get_lin_regress_params(initialPoints, 0, MIN_INIT_VALS, kx, bx, ky, by);
  526. cv::Rect_<track_t> rect0(kx * (MIN_INIT_VALS - 1) + bx, ky * (MIN_INIT_VALS - 1) + by, averageSize.x, averageSize.y);
  527. Point_t rectv0(kx, ky);
  528. switch (m_type)
  529. {
  530. case tracking::KalmanLinear:
  531. if (m_useAcceleration)
  532. CreateLinearAcceleration(rect0, rectv0);
  533. else
  534. CreateLinear(rect0, rectv0);
  535. break;
  536. }
  537. }
  538. }
  539. if (m_initialized)
  540. {
  541. cv::Mat measurement(4, 1, Mat_t(1));
  542. if (!dataCorrect)
  543. {
  544. measurement.at<track_t>(0) = m_lastRectResult.x; // update using prediction
  545. measurement.at<track_t>(1) = m_lastRectResult.y;
  546. measurement.at<track_t>(2) = m_lastRectResult.width;
  547. measurement.at<track_t>(3) = m_lastRectResult.height;
  548. }
  549. else
  550. {
  551. measurement.at<track_t>(0) = static_cast<track_t>(rect.x); // update using measurements
  552. measurement.at<track_t>(1) = static_cast<track_t>(rect.y);
  553. measurement.at<track_t>(2) = static_cast<track_t>(rect.width);
  554. measurement.at<track_t>(3) = static_cast<track_t>(rect.height);
  555. }
  556. // Correction
  557. cv::Mat estimated;
  558. switch (m_type)
  559. {
  560. case tracking::KalmanLinear:
  561. {
  562. estimated = m_linearKalman.correct(measurement);
  563. m_lastRectResult.x = estimated.at<track_t>(0); //update using measurements
  564. m_lastRectResult.y = estimated.at<track_t>(1);
  565. m_lastRectResult.width = estimated.at<track_t>(2);
  566. m_lastRectResult.height = estimated.at<track_t>(3);
  567. // Inertia correction
  568. if (!m_useAcceleration)
  569. {
  570. track_t currDist = sqrtf(sqr(estimated.at<track_t>(0) - rect.x) + sqr(estimated.at<track_t>(1) - rect.y) + sqr(estimated.at<track_t>(2) - rect.width) + sqr(estimated.at<track_t>(3) - rect.height));
  571. if (currDist > m_lastDist)
  572. {
  573. m_deltaTime = std::min(m_deltaTime + m_deltaStep, m_deltaTimeMax);
  574. }
  575. else
  576. {
  577. m_deltaTime = std::max(m_deltaTime - m_deltaStep, m_deltaTimeMin);
  578. }
  579. m_lastDist = currDist;
  580. m_linearKalman.transitionMatrix.at<track_t>(0, 4) = m_deltaTime;
  581. m_linearKalman.transitionMatrix.at<track_t>(1, 5) = m_deltaTime;
  582. m_linearKalman.transitionMatrix.at<track_t>(2, 6) = m_deltaTime;
  583. m_linearKalman.transitionMatrix.at<track_t>(3, 7) = m_deltaTime;
  584. }
  585. break;
  586. }
  587. }
  588. }
  589. else
  590. {
  591. if (dataCorrect)
  592. {
  593. m_lastRectResult.x = static_cast<track_t>(rect.x);
  594. m_lastRectResult.y = static_cast<track_t>(rect.y);
  595. m_lastRectResult.width = static_cast<track_t>(rect.width);
  596. m_lastRectResult.height = static_cast<track_t>(rect.height);
  597. }
  598. }
  599. return cv::Rect(static_cast<int>(m_lastRectResult.x), static_cast<int>(m_lastRectResult.y), static_cast<int>(m_lastRectResult.width), static_cast<int>(m_lastRectResult.height));
  600. }
  601. //---------------------------------------------------------------------------
  602. Rect3D TKalmanFilter::GetRect3DPrediction()
  603. {
  604. if (m_initialized)
  605. {
  606. cv::Mat prediction;
  607. switch (m_type)
  608. {
  609. case tracking::KalmanLinear:
  610. prediction = m_linearKalman.predict();
  611. break;
  612. }
  613. m_lastRect3DResult = Rect3D(cv::Point3f(prediction.at<track_t>(0), prediction.at<track_t>(1), prediction.at<track_t>(2)), Size3D(prediction.at<track_t>(3), prediction.at<track_t>(4), prediction.at<track_t>(5)), prediction.at<track_t>(6));
  614. }
  615. else
  616. {
  617. }
  618. return m_lastRect3DResult;
  619. }
  620. //---------------------------------------------------------------------------
  621. Rect3D TKalmanFilter::Update(Rect3D rect, bool dataCorrect)
  622. {
  623. if (!m_initialized)
  624. {
  625. if (m_initialRect3Ds.size() < MIN_INIT_VALS)
  626. {
  627. if (dataCorrect)
  628. {
  629. m_initialRect3Ds.push_back(rect);
  630. m_lastRect3DResult.center.x = static_cast<track_t>(rect.center.x);
  631. m_lastRect3DResult.center.y = static_cast<track_t>(rect.center.y);
  632. m_lastRect3DResult.center.z = static_cast<track_t>(rect.center.z);
  633. m_lastRect3DResult.size.width = static_cast<track_t>(rect.size.width);
  634. m_lastRect3DResult.size.height = static_cast<track_t>(rect.size.height);
  635. m_lastRect3DResult.size.length = static_cast<track_t>(rect.size.length);
  636. m_lastRect3DResult.yaw = static_cast<track_t>(rect.yaw);
  637. }
  638. }
  639. if (m_initialRect3Ds.size() == MIN_INIT_VALS)
  640. {
  641. std::vector<Point_t> initialPoints;
  642. cv::Point3f averageSize(0, 0, 0);
  643. float averageZ = 0;
  644. float averageYaw = 0;
  645. for (const auto& r : m_initialRect3Ds)
  646. {
  647. initialPoints.emplace_back(static_cast<track_t>(r.center.x), static_cast<track_t>(r.center.y));
  648. averageZ += r.center.z;
  649. averageSize.x += r.size.width;
  650. averageSize.y += r.size.height;
  651. averageSize.z += r.size.length;
  652. averageYaw += r.yaw;
  653. }
  654. averageZ /= MIN_INIT_VALS;
  655. averageSize.x /= MIN_INIT_VALS;
  656. averageSize.y /= MIN_INIT_VALS;
  657. averageSize.z /= MIN_INIT_VALS;
  658. averageYaw /= MIN_INIT_VALS;
  659. track_t kx = 0;
  660. track_t bx = 0;
  661. track_t ky = 0;
  662. track_t by = 0;
  663. get_lin_regress_params(initialPoints, 0, MIN_INIT_VALS, kx, bx, ky, by);
  664. Rect3D rect0(cv::Point3f(kx * (MIN_INIT_VALS - 1) + bx, ky * (MIN_INIT_VALS - 1) + by, averageZ), Size3D(averageSize.x, averageSize.y,averageSize.z), averageYaw);
  665. cv::Point3f rectv0(kx, ky, 0);
  666. switch (m_type)
  667. {
  668. case tracking::KalmanLinear:
  669. if (m_useAcceleration)
  670. CreateLinearAcceleration3D(rect0, rectv0);
  671. else
  672. CreateLinear3D(rect0, rectv0);
  673. break;
  674. }
  675. }
  676. }
  677. if (m_initialized)
  678. {
  679. cv::Mat measurement(7, 1, Mat_t(1));
  680. if (!dataCorrect)
  681. {
  682. measurement.at<track_t>(0) = m_lastRect3DResult.center.x; // update using prediction
  683. measurement.at<track_t>(1) = m_lastRect3DResult.center.y;
  684. measurement.at<track_t>(2) = m_lastRect3DResult.center.z;
  685. measurement.at<track_t>(3) = m_lastRect3DResult.size.width;
  686. measurement.at<track_t>(4) = m_lastRect3DResult.size.height;
  687. measurement.at<track_t>(5) = m_lastRect3DResult.size.length;
  688. measurement.at<track_t>(6) = m_lastRect3DResult.yaw;
  689. }
  690. else
  691. {
  692. measurement.at<track_t>(0) = static_cast<track_t>(rect.center.x); // update using measurements
  693. measurement.at<track_t>(1) = static_cast<track_t>(rect.center.y);
  694. measurement.at<track_t>(2) = static_cast<track_t>(rect.center.z);
  695. measurement.at<track_t>(3) = static_cast<track_t>(rect.size.width);
  696. measurement.at<track_t>(4) = static_cast<track_t>(rect.size.height);
  697. measurement.at<track_t>(5) = static_cast<track_t>(rect.size.length);
  698. measurement.at<track_t>(6) = static_cast<track_t>(rect.yaw);
  699. }
  700. // Correction
  701. cv::Mat estimated;
  702. switch (m_type)
  703. {
  704. case tracking::KalmanLinear:
  705. {
  706. estimated = m_linearKalman.correct(measurement);
  707. m_lastRect3DResult.center.x = estimated.at<track_t>(0); //update using measurements
  708. m_lastRect3DResult.center.y = estimated.at<track_t>(1);
  709. m_lastRect3DResult.center.z = estimated.at<track_t>(2);
  710. m_lastRect3DResult.size.width = estimated.at<track_t>(3);
  711. m_lastRect3DResult.size.height = estimated.at<track_t>(4);
  712. m_lastRect3DResult.size.length = estimated.at<track_t>(5);
  713. m_lastRect3DResult.yaw = estimated.at<track_t>(6);
  714. // Inertia correction
  715. if (!m_useAcceleration)
  716. {
  717. track_t currDist = sqrtf(sqr(estimated.at<track_t>(0) - rect.center.x) + sqr(estimated.at<track_t>(1) - rect.center.y)+ sqr(estimated.at<track_t>(2) - rect.center.z) + sqr(estimated.at<track_t>(3) - rect.size.width) + sqr(estimated.at<track_t>(4) - rect.size.height) + sqr(estimated.at<track_t>(5) - rect.size.length) + sqr(estimated.at<track_t>(6) - rect.yaw));
  718. if (currDist > m_lastDist)
  719. {
  720. m_deltaTime = std::min(m_deltaTime + m_deltaStep, m_deltaTimeMax);
  721. }
  722. else
  723. {
  724. m_deltaTime = std::max(m_deltaTime - m_deltaStep, m_deltaTimeMin);
  725. }
  726. m_lastDist = currDist;
  727. m_linearKalman.transitionMatrix.at<track_t>(0, 7) = m_deltaTime;
  728. m_linearKalman.transitionMatrix.at<track_t>(1, 8) = m_deltaTime;
  729. m_linearKalman.transitionMatrix.at<track_t>(2, 9) = m_deltaTime;
  730. m_linearKalman.transitionMatrix.at<track_t>(3, 10) = m_deltaTime;
  731. m_linearKalman.transitionMatrix.at<track_t>(4, 11) = m_deltaTime;
  732. m_linearKalman.transitionMatrix.at<track_t>(5, 12) = m_deltaTime;
  733. m_linearKalman.transitionMatrix.at<track_t>(6, 13) = m_deltaTime;
  734. }
  735. break;
  736. }
  737. }
  738. }
  739. else
  740. {
  741. if (dataCorrect)
  742. {
  743. m_lastRect3DResult.center.x = static_cast<track_t>(rect.center.x);
  744. m_lastRect3DResult.center.y = static_cast<track_t>(rect.center.y);
  745. m_lastRect3DResult.center.z = static_cast<track_t>(rect.center.z);
  746. m_lastRect3DResult.size.width = static_cast<track_t>(rect.size.width);
  747. m_lastRect3DResult.size.height = static_cast<track_t>(rect.size.height);
  748. m_lastRect3DResult.size.length = static_cast<track_t>(rect.size.length);
  749. m_lastRect3DResult.yaw = static_cast<track_t>(rect.yaw);
  750. }
  751. }
  752. return m_lastRect3DResult;
  753. }
  754. //---------------------------------------------------------------------------
  755. cv::Vec<track_t, 2> TKalmanFilter::GetVelocity() const
  756. {
  757. cv::Vec<track_t, 2> res(0, 0);
  758. if (m_initialized)
  759. {
  760. switch (m_type)
  761. {
  762. case tracking::KalmanLinear:
  763. {
  764. if (m_linearKalman.statePre.rows > 3)
  765. {
  766. int indX = 2;
  767. int indY = 3;
  768. if (m_linearKalman.statePre.rows > 5)
  769. {
  770. indX = 4;
  771. indY = 5;
  772. if (m_linearKalman.statePre.rows > 8)
  773. {
  774. indX = 7;
  775. indY = 8;
  776. }
  777. res[0] = m_linearKalman.statePre.at<track_t>(indX);
  778. res[1] = m_linearKalman.statePre.at<track_t>(indY);
  779. }
  780. }
  781. break;
  782. }
  783. }
  784. }
  785. return res;
  786. }