locate_algorithm.cpp 90 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830283128322833283428352836283728382839284028412842284328442845284628472848284928502851285228532854285528562857285828592860286128622863286428652866286728682869287028712872287328742875287628772878287928802881288228832884288528862887288828892890289128922893289428952896289728982899290029012902290329042905290629072908290929102911291229132914291529162917291829192920292129222923292429252926292729282929293029312932293329342935293629372938293929402941294229432944294529462947294829492950295129522953295429552956295729582959296029612962296329642965296629672968296929702971297229732974297529762977297829792980298129822983298429852986298729882989299029912992299329942995299629972998299930003001300230033004300530063007300830093010301130123013301430153016301730183019302030213022302330243025302630273028302930303031303230333034303530363037303830393040304130423043304430453046304730483049305030513052305330543055305630573058305930603061306230633064306530663067306830693070307130723073307430753076307730783079308030813082308330843085308630873088308930903091309230933094309530963097309830993100310131023103310431053106310731083109311031113112311331143115311631173118311931203121312231233124312531263127312831293130313131323133313431353136313731383139314031413142314331443145314631473148314931503151315231533154315531563157315831593160316131623163316431653166316731683169317031713172317331743175317631773178317931803181318231833184318531863187318831893190319131923193319431953196319731983199320032013202320332043205320632073208320932103211321232133214321532163217321832193220322132223223322432253226322732283229323032313232323332343235323632373238323932403241324232433244324532463247324832493250325132523253325432553256325732583259326032613262326332643265326632673268326932703271327232733274327532763277327832793280328132823283328432853286328732883289329032913292329332943295329632973298329933003301330233033304330533063307330833093310331133123313331433153316331733183319332033213322332333243325332633273328332933303331333233333334333533363337333833393340334133423343334433453346334733483349335033513352335333543355335633573358335933603361336233633364336533663367336833693370337133723373337433753376337733783379338033813382338333843385338633873388338933903391339233933394339533963397339833993400340134023403340434053406340734083409341034113412341334143415341634173418341934203421342234233424342534263427342834293430343134323433343434353436343734383439344034413442344334443445344634473448344934503451345234533454345534563457345834593460346134623463346434653466346734683469347034713472347334743475347634773478347934803481348234833484348534863487348834893490349134923493349434953496349734983499350035013502350335043505350635073508350935103511351235133514351535163517351835193520352135223523352435253526352735283529353035313532353335343535353635373538353935403541354235433544354535463547354835493550
  1. #include "stdafx.h"
  2. #include "locate_algorithm.h"
  3. #include <math.h>
  4. #include <algorithm>
  5. #include "../common/matrix/_Matrix.h"
  6. #include "../common/Functions/Functions.h"
  7. #include "ProcessRemodule.h"
  8. #include "log_process_module.h"
  9. using namespace std;
  10. // ant 天线, tar定位目标, d 测距
  11. void AOAformula( struct _coordinate* ant, double angle, struct _coordinate* tar, bool is_line /*= false*/, double z_offset /*= 0*/ )
  12. {
  13. double aa = sin(angle);
  14. double bb = cos(angle);
  15. double h = ant->z;
  16. if(is_line){
  17. h = z_offset;
  18. }
  19. tar->x = ant->x + sqrt(ant->d * ant->d - h * h) * cos(angle);
  20. tar->y = ant->y + sqrt(ant->d * ant->d - h * h) * sin(angle);
  21. tar->z = ant->z; // 高度
  22. tar->a = ant->a; // 角度
  23. tar->t = ant->t; // 时间戳
  24. tar->d_offset = ant->d_offset;
  25. tar->reader_id = ant->reader_id;
  26. tar->d = ant->d;
  27. }
  28. // ant 天线, last 上次坐标
  29. double Distanceformula( struct _coordinate* ant, struct _coordinate* last )
  30. {
  31. double x0,y0;
  32. double tan_diff = tan(last->a) - tan(ant->a);
  33. double ret = 0;
  34. if(tan_diff == 0){
  35. ret = sqrt((ant->x - last->x) * (ant->x - last->x) + (ant->y - last->y) * (ant->y - last->y));
  36. }else{
  37. x0 = (last->x * tan(last->a) - ant->x * tan(ant->a) + ant->y - last->y) / (tan(last->a) - tan(ant->a));
  38. y0 = (((last->x - ant->x) * tan(last->a) * tan(ant->a)) + ant->y * tan(last->a) - last->y * tan(ant->a)) / (tan(last->a) - tan(ant->a));
  39. ret = sqrt((x0 - ant->x) * (x0 - ant->x) + (y0 - ant->y) * (y0 - ant->y))
  40. + sqrt((x0 - last->x) * (x0 - last->x) + (y0 - last->y) * (y0 - last->y));
  41. }
  42. // n = Nodeformula( n1, n2);
  43. return ret;
  44. }
  45. // 接收到的新数据
  46. void algorithm_locate( struct _coordinate** coor_list, int coor_size, struct _coordinate* hist, struct _coordinate* dest, double dist_limit, double z_offset)
  47. {
  48. double d0; //历史定位坐标与实时定位坐标之间的距离;
  49. double dN; //利用速率V和时间戳理论计算出行进距离;
  50. double Angle1,Angle2,Angle3; //三个坐标时,任何2坐标的角度
  51. double Angle; //两个坐标之间的角度
  52. double dLong;
  53. double dLenth;
  54. int position; //定位点位置
  55. switch(coor_size){
  56. case 1: // 1条定位数据,5.3.1
  57. {
  58. if(coor_list[0]->d > dist_limit){ // 如果测距的距离大于某个固定值,此时判断此位置与天线角度一致
  59. //coor_list[0]->d -= coor_list[0]->d_offset;
  60. AOAformula(coor_list[0], coor_list[0]->a, dest, true, z_offset);
  61. }else { // 测距距离小于某值
  62. //利用速率V和时间戳理论计算出行进距离
  63. dN = coor_list[0]->v * (coor_list[0]->t - hist->t);
  64. //计算实时定位分站(天线)坐标和历史定位坐标之间的距离
  65. d0 = Distanceformula(coor_list[0], hist);
  66. if(hist->x - coor_list[0]->x == 0){
  67. Angle = (coor_list[0]->y >= hist->y) ? 90 : 270;
  68. Angle = Angle * M_PI / 180;
  69. }else{
  70. //计算历史坐标与实时分站坐标的角度,用来表明历史坐标位于定位坐标的方向
  71. Angle = atan((hist->y - coor_list[0]->y)/(hist->x - coor_list[0]->x));
  72. }
  73. // 与历史坐标同向
  74. //如果定位距离小于两个坐标之间的距离,说明定位位置在两个坐标之内
  75. if(d0 >= coor_list[0]->d){
  76. if(dN >= coor_list[0]->d){
  77. Angle += M_PI; //与历史坐标反方向
  78. }
  79. }else{
  80. //如果定位距离大于两个坐标之间的距离,说明定位位置在两个坐标之外,接下来只要判断在哪个yu坐标的位置即可
  81. //如果理论进行距离大于两者坐标之间的距离,则说明与历史坐标在同一方向,反之在反方向
  82. if(dN < coor_list[0]->d){
  83. Angle += M_PI; //与历史坐标反方向
  84. }
  85. }
  86. //coor_list[0]->d -= coor_list[0]->d_offset;
  87. AOAformula(coor_list[0], Angle, dest, true, z_offset); //最终定位坐标计算
  88. }
  89. break;
  90. }
  91. case 2: //2条定位数据,5.3.2
  92. {
  93. //2条定位数据中分站(天线)坐标与历史定位坐标之间的角度,以此推出这三个坐标是否在一条直线上或者呈三角形
  94. if(coor_list[0]->x == hist->x ){
  95. Angle1 = (coor_list[0]->y == hist->y) ? 0 : ((coor_list[0]->y > hist->y) ? 90 : 270);
  96. Angle1 = Angle1 * M_PI / 180;
  97. }else{
  98. Angle1 = atan((coor_list[0]->y - hist->y) / (coor_list[0]->x - hist->x));
  99. }
  100. if(coor_list[1]->x == hist->x ){
  101. Angle2 = (coor_list[1]->y == hist->y) ? 0 : ((coor_list[1]->y > hist->y)? 90 : 270);
  102. Angle2 = Angle2 * M_PI / 180;
  103. }else{
  104. Angle2 = atan((coor_list[1]->y - hist->y) / (coor_list[1]->x - hist->x));
  105. }
  106. if(coor_list[1]->x == coor_list[0]->x ){
  107. Angle3 = (coor_list[1]->y >= coor_list[0]->y) ? 90 : 270; // 必须在不同点
  108. Angle3 = Angle3 * M_PI / 180;
  109. }else{
  110. Angle3 = atan((coor_list[1]->y - coor_list[0]->y) / (coor_list[1]->x - coor_list[0]->x));
  111. if(Angle3 > M_PI){
  112. Angle3 -= M_PI;
  113. }
  114. /*if(coor_list[0]->x > coor_list[1]->x)
  115. Angle3 += M_PI;*/
  116. }
  117. //2个角度相同或者相差180度说明在一条直线上
  118. //if(coor_list[0]->a == coor_list[1]->a || coor_list[0]->a == coor_list[1]->a + M_PI || coor_list[0]->a == coor_list[1]->a - M_PI){
  119. if(Angle1 == Angle2 || Angle1 == Angle2 + M_PI || Angle1 == Angle2 - M_PI){
  120. //取两个距离之间的最大值,保存最短距离的数据作为定位点
  121. dLong = max(coor_list[0]->d, coor_list[1]->d);
  122. position = (coor_list[0]->d > coor_list[1]->d) ? 1 : 0;
  123. dLenth = Distanceformula(coor_list[0], coor_list[1]);
  124. if(dLong > dLenth){
  125. Angle = Angle3 + ((position + 1)%2) * M_PI;
  126. }else{
  127. Angle = Angle3 + position * M_PI;
  128. }
  129. //2个坐标的距离小于其2个距离中的最大值,则说明定位点在2个坐标的外侧,
  130. //且在距离短的那条作为定位坐标,定位角度与此相同
  131. //2个坐标的距离大于这两条距离中的任何一个,则说明定位点在两个坐标之间,取任意一条作为定位点,定位角度与其角度相反(+180)
  132. // AOAformula(coor_list[position], Angle, dest); //最终定位坐标计算
  133. // 根据测量的结果,调整距离d
  134. // x2 - (d2 - ((d1+d2) - (x2-x1))/2) 靠近x2 x2 > x1,d1 > d2
  135. // x1 + (d1 - ((d1+d2) - (x2-x1))/2) 靠近x1 x2 > x1,d1 < d2
  136. // 化简得到 (x1 + x2 + d1 - d2)/2
  137. // 显示的d = ((d1 + d2) - abs(x2-x1))/2
  138. double d_offset;
  139. if(fabs(coor_list[0]->x - coor_list[1]->x) > fabs(coor_list[0]->y - coor_list[1]->y)){ // 水平方向
  140. d_offset = (coor_list[0]->d + coor_list[1]->d - fabs(coor_list[0]->x - coor_list[1]->x))/2;
  141. }else{ // 垂直方向
  142. d_offset = (coor_list[0]->d + coor_list[1]->d - fabs(coor_list[0]->y - coor_list[1]->y))/2;
  143. }
  144. dest->d_offset = d_offset;
  145. coor_list[0]->d_offset = d_offset;
  146. coor_list[1]->d_offset = d_offset;
  147. coor_list[position]->d -= d_offset;
  148. AOAformula(coor_list[position], Angle, dest, true, z_offset);
  149. }
  150. //岔路口形式
  151. else {
  152. dLenth = Distanceformula(coor_list[0], coor_list[1]); //2个天线坐标之间的距离
  153. if(dLenth > (coor_list[0]->d + coor_list[1]->d)) { //取一条作为定位参考数据,采用1条定位数据算法,参考5.3.1
  154. //如果测距的距离大于某个固定值,此时判断此位置与天线角度一致,
  155. if(coor_list[0]->d > dist_limit){
  156. AOAformula(coor_list[0], coor_list[0]->a, dest); //最终定位坐标计算
  157. }
  158. //测距距离小于某值
  159. else{
  160. dN = coor_list[0]->v * ( hist->t - coor_list[0]->t); //利用速率V和时间戳理论计算出行进距离;
  161. d0 = Distanceformula(coor_list[0], hist); //计算实时定位分站(天线)坐标和历史定位坐标之间的距离
  162. if(hist->x == coor_list[0]->x){
  163. Angle = (hist->y >= coor_list[0]->y) ? 90 : 270; //计算历史坐标与实时分站坐标的角度,用来表明历史坐标位于定位坐标的方向
  164. Angle = Angle * M_PI / 180;
  165. }else{
  166. Angle = atan((hist->y - coor_list[0]->y) / (hist->x - coor_list[0]->x));
  167. }
  168. //如果定位距离小于两个坐标之间的距离,说明定位位置在两个坐标之内
  169. if(d0 > coor_list[0]->d){
  170. if(dN >= coor_list[0]->d){
  171. Angle += M_PI; //与历史坐标反方向
  172. }
  173. }
  174. else{
  175. //如果定位距离大于两个坐标之间的距离,说明定位位置在两个坐标之外,接下来只要判断在哪个坐标的位置即可
  176. //如果理论进行距离大于两者坐标之间的距离,则说明与历史坐标在同一方向,反之在反方向
  177. if(dN < coor_list[0]->d){
  178. Angle += M_PI;
  179. }
  180. }
  181. AOAformula(coor_list[0], Angle, dest, true, z_offset); //最终定位坐标计算
  182. }
  183. }
  184. else { //采用2条定位数据算法,采用5.3.2
  185. //取两个距离之间的最大值,保存最短距离的数据作为定位点
  186. dLong = max(coor_list[0]->d, coor_list[1]->d);
  187. position = (coor_list[0]->d > coor_list[1]->d) ? 1 : 0;
  188. dLenth = Distanceformula(coor_list[0], coor_list[1]);
  189. //2个坐标的距离小于其2个距离中的最大值,则说明定位点在2个坐标的外侧,
  190. //且在距离短的那条作为定位坐标,定位角度与此相同
  191. //2个坐标的距离大于这两条距离中的任何一个,则说明定位点在两个坐标之间,取任意一条作为定位点,定位角度与其角度相反
  192. if(dLong > dLenth){
  193. Angle = Angle3 + ((position + 1)%2) * M_PI;
  194. }else{
  195. Angle = Angle3 + position * M_PI;
  196. }
  197. AOAformula(coor_list[position], Angle, dest, true, z_offset); //最终定位坐标计算
  198. }
  199. }
  200. break;
  201. }
  202. case 3: // 3条定位数据,三点定位法
  203. {
  204. // 在同一直线上不能使用三点定位
  205. locatebycordinate(coor_list, dest);
  206. break;
  207. }
  208. case 4: // 4条数据,取前3条
  209. {
  210. locatebycordinate(coor_list, dest);
  211. break;
  212. }
  213. default:
  214. {
  215. dest->x = -1;
  216. dest->y = -1;
  217. dest->a = 0;
  218. break;
  219. }
  220. }
  221. }
  222. void locatebycordinate( struct _coordinate** coor_list, struct _coordinate* dest )
  223. {
  224. double h = 0;
  225. int i = 0;
  226. int j = 0;
  227. double temp = 0;
  228. double x[3]; // 坐标
  229. double y[3]; // 坐标
  230. double z[3];
  231. double d[3];
  232. double p[2];
  233. _Matrix pi;
  234. matrix_set(&pi,3,2);
  235. matrix_init(&pi);
  236. for(i=0;i<3;i++){
  237. x[i] = coor_list[i]->x;
  238. y[i] = coor_list[i]->y;
  239. z[i] = coor_list[i]->z;
  240. d[i] = coor_list[i]->d; // d * CM_PIX * 100;
  241. }
  242. //获得坐标
  243. matrix_write(&pi,0,0,x[0]);
  244. matrix_write(&pi,0,1,y[0]);
  245. h = z[0];
  246. d[0] = sqrt(d[0] * d[0] - h * h);
  247. matrix_write(&pi,1,0,x[1]);
  248. matrix_write(&pi,1,1,y[1]);
  249. h = z[1];
  250. d[1] = sqrt(d[1] * d[1]- h * h);
  251. matrix_write(&pi,2,0,x[2]);
  252. matrix_write(&pi,2,1,y[2]);
  253. h = z[2];
  254. d[2] = sqrt(d[2] * d[2] - h * h);
  255. if (locate_func(&pi,d,p) > 0)
  256. {
  257. dest->x = p[0];
  258. dest->y = p[1];
  259. }else
  260. {
  261. dest->x = -1;
  262. dest->y = -1;
  263. }
  264. }
  265. void algorithm_locate_ex( _coordinate** coor_list, int coor_size, _coordinate* hist, _coordinate* dest, double dist_limit )
  266. {
  267. switch (coor_size)
  268. {
  269. case 1:
  270. {
  271. if(!hist){ // 初次计算,只有一个分站数据的时候跳出,不计算坐标
  272. return;
  273. }
  274. // 获取历史节点与当前分站的角度,判断是否仍在同一直线
  275. double angle = getAngleByAngle(coor_list[0], hist);
  276. if(0 == tan(angle)){ // 在同一直线
  277. double dist_temp = getDistance(hist, coor_list[0]); // 获得上一次的距离
  278. // 获得移动的距离
  279. double dist_mv = fabs(dist_temp - coor_list[0]->d);
  280. // 沿着原移动方向运动
  281. calcCoordinate(hist, dist_mv, hist->a, dest); //
  282. }else{ // 有角度,找交点
  283. _coordinate* pt_cross = getCross(hist, coor_list[0]);
  284. double dist_temp = getDistance(pt_cross, coor_list[0]);
  285. if(dist_temp >= coor_list[0]->d){
  286. // 避免动画穿墙,先将该点定位到交点上
  287. dest->x = pt_cross->x;
  288. dest->y = pt_cross->y;
  289. dest->z = coor_list[0]->z;
  290. dest->a = getAngleByCoordinate(pt_cross, coor_list[0]); // 运行方向
  291. }else{ // 在当前方向上运行, 理论上不存在
  292. dist_temp = (getDistance(hist, pt_cross) - (coor_list[0]->d - dist_temp)); // 位移
  293. calcCoordinate(hist, dist_temp, hist->a, dest);
  294. }
  295. }
  296. break;
  297. }
  298. case 2:
  299. {
  300. double angle_anc = getAngleByAngle(coor_list[0], coor_list[1]);
  301. if(0 == tan(angle_anc)){ // 在同一直线
  302. double dist_anc = getDistance(coor_list[0], coor_list[1]);
  303. double dist_offset;
  304. if(dist_anc < coor_list[0]->d || dist_anc < coor_list[1]->d ) {//分站到标识卡的距离大于分站之间的距离,同一侧
  305. dist_offset = fabs(coor_list[0]->d - coor_list[1]->d) - dist_anc;
  306. double angle_anc_run = getAngleByCoordinate(coor_list[0], coor_list[1]); // 分站2一侧
  307. if(coor_list[0]->d < coor_list[1]->d){ // 分站1一侧
  308. angle_anc_run += M_PI;
  309. }
  310. calcCoordinate(coor_list[0], coor_list[0]->d - dist_offset / 2, angle_anc_run, dest);
  311. }else{ //在分站之间
  312. dist_offset = coor_list[0]->d + coor_list[1]->d - dist_anc;
  313. calcCoordinate(coor_list[0], coor_list[0]->d + dist_offset / 2, getAngleByCoordinate(coor_list[0], coor_list[1]), dest);
  314. }
  315. }else{ // 有夹角, 一定在两分站之间
  316. if(hist){
  317. double angle_anc1 = getAngleByAngle(coor_list[0], hist);
  318. double angle_anc2 = getAngleByAngle(coor_list[1], hist);
  319. // 找到交点
  320. _coordinate* pt_cross = getCross(coor_list[0], coor_list[1]);
  321. // 交点到分站的距离
  322. double dist_anc1 = getDistance(pt_cross, coor_list[0]);
  323. double dist_anc2 = getDistance(pt_cross, coor_list[1]);
  324. if(0 == tan(angle_anc1)){
  325. if(dist_anc1 >= coor_list[0]->d){ // 以直线上为准
  326. calcCoordinate(hist, fabs(dist_anc1 - coor_list[0]->d), hist->a, dest);
  327. //calcCoordinate(coor_list[0], coor_list[0]->d, getAngleByCoordinate(coor_list[0], hist), dest);
  328. }else{ // 超过直接定位到交点上
  329. dest->x = pt_cross->x;
  330. dest->y = pt_cross->y;
  331. dest->a = getAngleByCoordinate(pt_cross, coor_list[1]);
  332. }
  333. }else if(0 == tan(angle_anc2)){ // 与分站1在同一直线
  334. if(dist_anc2 >= coor_list[1]->d){ // 以直线上为准
  335. calcCoordinate(hist, fabs(dist_anc2 - coor_list[1]->d), hist->a, dest);
  336. //calcCoordinate(coor_list[1], coor_list[1]->d, getAngleByCoordinate(coor_list[1], hist), dest);
  337. }else{ // 超过直接定位到交点上
  338. dest->x = pt_cross->x;
  339. dest->y = pt_cross->y;
  340. dest->a = getAngleByCoordinate(pt_cross, coor_list[0]);
  341. }
  342. }else { // 三叉口
  343. // 找到交点
  344. _coordinate* pt_cross1 = getCross(hist, coor_list[0]);
  345. _coordinate* pt_cross2 = getCross(hist, coor_list[1]);
  346. if(coor_list[0]->d >= coor_list[1]->d){
  347. dest->x = pt_cross2->x;
  348. dest->y = pt_cross2->y;
  349. dest->z = coor_list[1]->z;
  350. dest->a = coor_list[1]->a;
  351. }else{
  352. dest->x = pt_cross1->x;
  353. dest->y = pt_cross1->y;
  354. dest->z = coor_list[0]->z;
  355. dest->a = coor_list[0]->a;
  356. }
  357. }
  358. }else{ // 先定位到交点,避免穿墙
  359. _coordinate* pt_cross = getCross(coor_list[0], coor_list[1]);
  360. dest->x = pt_cross->x;
  361. dest->y = pt_cross->y;
  362. // 比较距离,定位到距离近的方向上
  363. if(fabs(getDistance(pt_cross, coor_list[0]) - coor_list[0]->d) >= fabs(getDistance(pt_cross, coor_list[1]) - coor_list[1]->d)){
  364. dest->z = coor_list[1]->z;
  365. dest->a = coor_list[1]->a;
  366. }else{
  367. dest->z = coor_list[0]->z;
  368. dest->a = coor_list[0]->a;
  369. }
  370. }
  371. }
  372. break;
  373. }
  374. case 3:
  375. {
  376. // 判断分站是否在同一直线上
  377. if(is_all_inline(coor_list, coor_size)){
  378. algorithm_locate_ex(coor_list, 2, hist, dest, dist_limit);
  379. break;
  380. }
  381. locatebycordinate(coor_list, dest);
  382. break;
  383. }
  384. case 4:
  385. {
  386. // 判断分站是否在同一直线上
  387. if(is_all_inline(coor_list, coor_size)){
  388. algorithm_locate_ex(coor_list, 2, hist, dest, dist_limit);
  389. break;
  390. }
  391. locatebycordinate(coor_list, dest);
  392. break;
  393. }
  394. default:
  395. break;
  396. }
  397. }
  398. double getAngleByAngle( _coordinate* p1, _coordinate* p2 )
  399. {
  400. if(p1->x == p2->x) return 0;
  401. if(p1->y == p2->y) return 0;
  402. return double(p2->a - p1->a);
  403. }
  404. double getAngleByCoordinate( _coordinate* p1, _coordinate* p2 )
  405. {
  406. if(p1->x == p2->x && p1->y == p2->y){ // 如果坐标相等,则跨越,不考虑折返情况
  407. return p1->a;
  408. }
  409. if(p1->x == p2->x){
  410. if(p1->y >= p2->y){
  411. return M_PI_2 * 3;
  412. }
  413. return M_PI_2;
  414. }
  415. if(p1->y == p2->y){
  416. if(p1->x > p2->x){
  417. return M_PI;
  418. }
  419. return 0;
  420. }
  421. return atan((p2->y - p1->y) / (p2->x - p1->x));
  422. }
  423. _coordinate* getCross( _coordinate* p1, _coordinate* p2 )
  424. {
  425. _coordinate* p = new _coordinate;
  426. double b1 = p1->y - tan(p1->a) * p1->x;
  427. double b2 = p2->y - tan(p2->a) * p2->x;
  428. if(p1->a == 0 || p1->a == M_PI){ // 第一个点水平
  429. p->y = p1->y;
  430. if(p2->a == M_PI_2 || p2->a == M_PI_2 * 3){ // 第二个点垂直
  431. p->x = p2->x;
  432. }else{ // 第二个点有角度
  433. p->x = (p->y - b2) / tan(p2->a);
  434. }
  435. return p;
  436. }else if(p1->a == M_PI_2 || p1->a == M_PI_2 * 3){ // 第一个点垂直
  437. p->x = p1->x;
  438. if(p2->a == 0 || p2->a == M_PI){ // 第二个点水平
  439. p->y = p2->y;
  440. }else{ // 第二个点有角度
  441. p->y = tan(p2->a) * p->x + b2;
  442. }
  443. return p;
  444. }else { // 有角度
  445. if(p2->a == 0 || p2->a == M_PI){ // 第二点水平
  446. p->y = p2->y;
  447. p->x = (p->y - b1) / tan(p1->a);
  448. return p;
  449. }else if(p2->a == M_PI_2 || p2->a == M_PI_2 * 3){ // 第二点垂直
  450. p->x = p2->x;
  451. p->y = p->x * tan(p1->a) + b1;
  452. return p;
  453. }
  454. }
  455. p->x = (b1-b2)/(tan(p2->a) - tan(p1->a));
  456. p->y = tan(p1->a) * p->x + b1;
  457. return p;
  458. }
  459. // 不考虑折返
  460. void calcCoordinate( _coordinate* p1, double dist, double angle, _coordinate* dest )
  461. {
  462. dest->x = p1->x;
  463. dest->y = p1->y;
  464. dest->z = p1->z;
  465. dest->x += dist * cos(angle);
  466. dest->y += dist * sin(angle);
  467. dest->a = angle;
  468. }
  469. double getDistance( _coordinate* p1, _coordinate *p2 )
  470. {
  471. return sqrt(pow(p1->x - p2->x, 2) + pow(p1->y - p2->y, 2));
  472. }
  473. bool is_all_inline( _coordinate** coor_list, int cnt )
  474. {
  475. bool ret = false;
  476. if(cnt < 3) return ret;
  477. if(coor_list[0]->x == coor_list[1]->x && coor_list[0]->x == coor_list[2]->x){
  478. ret = true;
  479. }else if(coor_list[0]->y == coor_list[1]->y && coor_list[0]->y == coor_list[2]->y){
  480. ret = true;
  481. }else{
  482. ret = false;
  483. }
  484. return ret;
  485. }
  486. POS* LocateAlgorithm::Pos(ReaderPathMap rpm,int sta_num,int ant,double dist,INFO_PRE info_pre)
  487. {
  488. POS* p = new POS;
  489. p->posx = 0;
  490. p->posy = 0;
  491. p->pos_radius = 999999.9; //inf在matlab中表示负无穷
  492. double d[2];
  493. double xcross[2];
  494. double ycross[2];
  495. int seg_num = 0;
  496. if(rpm.find(sta_num) == rpm.end())
  497. {
  498. return NULL;
  499. }
  500. seg_num = rpm.find(sta_num)->second->nRealCalcPoints - 1;
  501. if(seg_num < 0 ){
  502. return NULL;
  503. }
  504. int j = 0;
  505. for(int i =0;i<seg_num;i++){
  506. if(rpm.find(sta_num)->second->px[i]==-1||rpm.find(sta_num)->second->px[i+1]==-1)
  507. {
  508. continue;
  509. }
  510. SOLUTION *r = NULL;
  511. r = GetPos(rpm,sta_num,ant,dist,i);
  512. if(r == NULL){
  513. continue;
  514. }
  515. //x的列数count取2,是因为解最多为2
  516. int count = 2;
  517. for(j =0;j<2;j++){
  518. xcross[j] = r->x[j];
  519. ycross[j] = r->y[j];
  520. }
  521. for(j = 0;j < count;j++){
  522. double d1 = xcross[j] - rpm.find(info_pre.sta_num)->second->x[info_pre.ant];
  523. double d2 = ycross[j] - rpm.find(info_pre.sta_num)->second->y[info_pre.ant];
  524. d[j] = sqrt(pow(d1,2) + pow(d2,2));
  525. }
  526. for(j = 0;j < count;j++){
  527. d[j] = fabs(d[j] - info_pre.dist);
  528. }
  529. double c = 999999.9;
  530. int k = -1;
  531. for(j = 0;j < count;j++){
  532. if(d[j] < c){
  533. c = d[j];
  534. k = j;
  535. }
  536. }
  537. if(c<p->pos_radius){
  538. p->posx = xcross[k];
  539. p->posy = ycross[k];
  540. p->pos_radius = c;
  541. }
  542. }
  543. return p;
  544. }
  545. /*
  546. * TOF定位:遍历路段集和指定分站天线为圆心的坐标
  547. *
  548. * param
  549. * trpm tof地图路径集
  550. * reader_id 分站id
  551. * antenna_id 分站天线索引
  552. * dist 距离此分站天线的距离
  553. * refer_data 参考数据信息
  554. *
  555. * return
  556. * 返回定位的坐标,如果无返回nullptr
  557. *
  558. */
  559. std::shared_ptr<POS> LocateAlgorithm::tof_locate_1d(std::shared_ptr<TOFReaderPathMap> trpm,int reader_id,int antenna_idx,double dist,TOF_REFER_DATA refer_data)
  560. {
  561. std::shared_ptr<POS> p = std::make_shared<POS>();
  562. double d[MAX_READER_TDOA_PATH_NUMS] = {0};
  563. double xcross[MAX_READER_TDOA_PATH_NUMS] = {0};
  564. double ycross[MAX_READER_TDOA_PATH_NUMS] = {0};
  565. int seg_num = 0;
  566. if(trpm->find(reader_id) == trpm->end() || trpm->find(refer_data.nReaderId) == trpm->end())
  567. {
  568. return nullptr;
  569. }
  570. seg_num = trpm->find(reader_id)->second->nRealCalcPoints - 1;
  571. //判断路径集的点数是否大于两个点,两点才可以构成线段路径
  572. if(seg_num < 0 ){
  573. return nullptr;
  574. }
  575. int j = 0;
  576. int m = 0;
  577. int nCounts = 0;
  578. //对每一个线段所在的直线和以天线为圆心的圆求解,最多为两个解
  579. for(int i = 0;i < seg_num;i++){
  580. //如果线段两端点的坐标不存在,则寻找下一个线段
  581. if(trpm->find(reader_id)->second->px[i] == -1 || trpm->find(reader_id)->second->px[i+1] == -1)
  582. {
  583. continue;
  584. }
  585. //求线段所在直线和分站readerid的antenna_idx天线为圆心,半径为dist的圆的两个交点
  586. std::shared_ptr<SOLUTION> r = GetPos(trpm,reader_id,antenna_idx,dist,i);
  587. //如果解不存在,则继续寻找下一个线段
  588. if(r == nullptr || r->nCount <= 0){
  589. continue;
  590. }
  591. int k = -1;
  592. double c = 999999.9;
  593. //如果解存在,则保存两个解
  594. for(j = 0;j < r->nCount;j++){
  595. //保存解
  596. xcross[j] = r->x[j];
  597. ycross[j] = r->y[j];
  598. //计算解与参考分站中的天线坐标的距离distance
  599. double d1 = xcross[j] - trpm->find(refer_data.nReaderId)->second->x[refer_data.nAntennaIndex];
  600. double d2 = ycross[j] - trpm->find(refer_data.nReaderId)->second->y[refer_data.nAntennaIndex];
  601. double distance = sqrt(pow(d1,2) + pow(d2,2));
  602. d[j] = distance;
  603. }
  604. for (j = 0;j < r->nCount; j++)
  605. {
  606. if (d[j] < c)
  607. {
  608. c = d[j];
  609. k = j;
  610. }
  611. }
  612. if(c < p->pos_radius){
  613. p->posx = xcross[k];
  614. p->posy = ycross[k];
  615. p->pos_radius = c;
  616. }
  617. }
  618. return p;
  619. }
  620. std::shared_ptr<POS> LocateAlgorithm::TdoaLocate2d(std::shared_ptr<ReceiveDataMap> pRdm)
  621. {
  622. if (pRdm->size() < 3)
  623. {
  624. return nullptr;
  625. }
  626. double BS[4][2] = {0};
  627. double h[3] = {0};
  628. unsigned int i = 0;
  629. ReceiveDataMap::iterator first = pRdm->begin();
  630. for (ReceiveDataMap::iterator it = pRdm->begin();it != pRdm->end();++it)
  631. {
  632. BS[i][0] = it->second->x;
  633. BS[i][1] = it->second->y;
  634. if (i > 0)
  635. {
  636. h[i-1] = CFunctions::getDistance(it->first - first->first,CFunctions::TDOA);
  637. }
  638. i++;
  639. if (i > pRdm->size())
  640. {
  641. break;
  642. }
  643. }
  644. std::shared_ptr<POS> p = std::make_shared<POS>();
  645. double a = (BS[0][0]+BS[1][0]+BS[2][0])/2;
  646. double c1 = 1.49445;
  647. double c2 = 1.49445;
  648. int maxgen = 1000;
  649. unsigned int sizepop = 100;
  650. int Vmax = 100;
  651. int Vmin = -100;
  652. int popmax = 50000;
  653. int popmin = -50000;
  654. /*初始化参数*/
  655. /*产生初始粒子和速度*/
  656. unsigned int j = 0;
  657. double fitness[1][100] = {0};
  658. double qq1 = 0,qq2 = 0,re = 0;
  659. double pop[100][2] = {0},V[100][2] = {0};//sizepop = 100;
  660. for(i=0;i<sizepop;i++)
  661. {
  662. for(j = 0;j<2;j++)
  663. {
  664. pop[i][j] =a* (rand()+1)/(RAND_MAX+2.0) ;
  665. V[i][j] = 10*(rand()+1)/(RAND_MAX+2.0);
  666. }
  667. }
  668. for(i = 0;i < sizepop;i++)
  669. {
  670. qq1 = pop[i][0];
  671. qq2 = pop[i][1];
  672. re = pow((sqrt((qq1-BS[1][0])*(qq1-BS[1][0])+( qq2-BS[1][1])*( qq2-BS[1][1]))-sqrt((qq1-BS[0][0])*(qq1-BS[0][0])+( qq2-BS[0][1])*( qq2-BS[0][1]))-h[0]),2)
  673. +pow(( sqrt((qq1-BS[2][0])*(qq1-BS[2][0])+( qq2-BS[2][1])*( qq2-BS[2][1]))-sqrt((qq1-BS[0][0])*(qq1-BS[0][0])+( qq2-BS[0][1])*( qq2-BS[0][1]))-h[1]),2);
  674. fitness[0][i]=re;
  675. }
  676. /*产生初始粒子和速度*/
  677. /*个体极值和群体极值*/
  678. double t = 0;
  679. double fitnesscopy[1][100] = {0};
  680. for(i=0;i<sizepop;i++)
  681. {
  682. fitnesscopy[0][i] = fitness[0][i];
  683. }
  684. for(i=0;i<sizepop;i++)
  685. {
  686. for(j=i+1;j<sizepop;j++)
  687. {
  688. if (fitnesscopy[0][i]>fitnesscopy[0][j])
  689. {
  690. t = fitnesscopy[0][i];
  691. fitnesscopy[0][i]=fitnesscopy[0][j];
  692. fitnesscopy[0][j] = t;
  693. }
  694. }
  695. }//正确
  696. double bestfitness = fitnesscopy[0][0];
  697. int bestindex = 0,kk = 0;
  698. for(i=0;i<sizepop;i++)
  699. {
  700. if(fitness[0][i]==bestfitness)
  701. {
  702. kk = i;
  703. }
  704. }
  705. bestindex = kk;//正确
  706. double zbest[1][2] = {0};
  707. for (j = 0;j<2;j++)
  708. {
  709. zbest[0][j] = pop[bestindex][j];//全局最佳
  710. }
  711. double gbest[100][2] = {0};
  712. for(i = 0;i<100;i++)
  713. {
  714. for (j = 0;j<2;j++)
  715. {
  716. gbest[i][j] = pop[i][j];//个体最佳
  717. }
  718. }
  719. double fitnessgbest[1][100] = {0};
  720. double fitnesszbest;
  721. for(i = 0;i<100;i++)
  722. {
  723. fitnessgbest[0][i]=fitness[0][i];
  724. }
  725. fitnesszbest = bestfitness;
  726. /*个体极值和群体极值*/
  727. /*迭代寻优*/
  728. double w;
  729. double ws = 0.9;
  730. double we = 0.4;
  731. for(i = 0;i<maxgen;i++)
  732. {
  733. w = ws - (ws-we)*(i/maxgen)*(i/maxgen);
  734. for(j=0;j<sizepop;j++)
  735. {
  736. V[j][0]= w*V[j][0]+c1*(rand()+1)/(RAND_MAX+2.0)*(gbest[j][0]-pop[j][0])+c2*(rand()+1)/(RAND_MAX+2.0)*(zbest[j][0]-pop[j][0]);
  737. V[j][1]= w*V[j][1]+c1*(rand()+1)/(RAND_MAX+2.0)*(gbest[j][1]-pop[j][1])+c2*(rand()+1)/(RAND_MAX+2.0)*(zbest[j][1]-pop[j][1]);
  738. if(V[j][0]>Vmax)
  739. {
  740. V[j][0]=Vmax;
  741. }
  742. if(V[j][1]>Vmax)
  743. {
  744. V[j][1]=Vmax;
  745. }
  746. if(V[j][0]<Vmin)
  747. {
  748. V[j][0]=Vmin;
  749. }
  750. if(V[j][1]<Vmin)
  751. {
  752. V[j][1]=Vmin;
  753. }
  754. //种群更新//
  755. pop[j][0]=pop[j][0]+ V[j][0];
  756. pop[j][1]= pop[j][1]+V[j][1];
  757. if(pop[j][0]>popmax)
  758. {
  759. pop[j][0]=popmax;
  760. }
  761. if(pop[j][1]>popmax)
  762. {
  763. pop[j][1]=popmax;
  764. }
  765. if(pop[j][0]<popmin)
  766. {
  767. pop[j][0]=popmin;
  768. }
  769. if(pop[j][1]<popmin)
  770. {
  771. pop[j][1]=popmin;
  772. }
  773. //种群更新//
  774. double qq3,qq4,re1;
  775. qq3 = pop[j][0];
  776. qq4 = pop[j][1];
  777. re1=pow((sqrt((qq3-BS[1][0])*(qq3-BS[1][0])+( qq4-BS[1][1])*( qq4-BS[1][1]))-sqrt((qq3-BS[0][0])*(qq3-BS[0][0])+( qq4-BS[0][1])*( qq4-BS[0][1]))-h[0]),2)
  778. +pow(( sqrt((qq3-BS[2][0])*(qq3-BS[2][0])+( qq4-BS[2][1])*( qq4-BS[2][1]))-sqrt((qq3-BS[0][0])*(qq3-BS[0][0])+( qq4-BS[0][1])*( qq4-BS[0][1]))-h[1]),2);
  779. fitness[0][j]=re1;
  780. }//正确
  781. for(j=0;j<sizepop;j++)
  782. {
  783. if( fitness[0][j] < fitnessgbest[0][j])
  784. {
  785. gbest[j][0] = pop[j][0];
  786. gbest[j][1] = pop[j][1];
  787. fitnessgbest[0][j] = fitness[0][j];
  788. }//个体最优更新
  789. if(fitness[0][j] < fitnesszbest)
  790. {
  791. zbest[0][0] = pop[j][0];
  792. zbest[0][1] = pop[j][1];
  793. fitnesszbest = fitness[0][j];
  794. }//群体最优更新
  795. }
  796. }
  797. p->posx = zbest[0][0];
  798. p->posy = zbest[0][1];
  799. return p;
  800. }
  801. std::shared_ptr<POS> LocateAlgorithm::TdoaLocate3d(std::shared_ptr<ReceiveDataMap> pRdm)
  802. {
  803. if (pRdm->size() < 3)
  804. {
  805. return nullptr;
  806. }
  807. std::shared_ptr<POS> p = std::make_shared<POS>();
  808. double BS[4][3] = {0};
  809. double h[3] = {0};
  810. int i = 0;
  811. ReceiveDataMap::iterator first = pRdm->begin();
  812. for (ReceiveDataMap::iterator it = pRdm->begin();it != pRdm->end();++it)
  813. {
  814. BS[i][0] = it->second->x;
  815. BS[i][1] = it->second->y;
  816. BS[i][2] = it->second->z;
  817. if (i > 0)
  818. {
  819. h[i-1] = CFunctions::getDistance(it->first - first->first,CFunctions::TDOA);
  820. }
  821. i++;
  822. if (i > pRdm->size())
  823. {
  824. break;
  825. }
  826. }
  827. /*初始化参数*/
  828. double a = (BS[0][0]+BS[1][0]+BS[2][0]+BS[3][0])/2;
  829. double c1=1.49445;
  830. double c2=1.49445;
  831. int maxgen = 1000;
  832. int sizepop = 100;
  833. int Vmax = 10;
  834. int Vmin = -10;
  835. int popmax = 50;
  836. int popmin = -50;
  837. /*初始化参数*/
  838. /*产生初始粒子和速度*/
  839. int j;
  840. double fitness[1][100] = {0};
  841. double qq1,qq2,qq5,re;
  842. double pop[100][3] = {0},V[100][3] = {0};//sizepop = 100;
  843. //srand((unsigned)time(NULL));
  844. for(i=0;i<sizepop;i++)
  845. {
  846. for(j = 0;j<3;j++)
  847. {
  848. //pop[i][j] =a* (rand()+1)/(RAND_MAX+2.0) ;
  849. //V[i][j] = 10*(rand()+1)/(RAND_MAX+2.0);
  850. pop[i][j] =a* (( rand())/double(RAND_MAX));
  851. V[i][j] = 10* ( (rand())/double(RAND_MAX));
  852. }
  853. }//正确
  854. for(i=0;i<sizepop;i++)
  855. {
  856. qq1= pop[i][0];
  857. qq2= pop[i][1];
  858. qq5 = pop[i][2];
  859. re=pow(( sqrt((qq1-BS[1][0])*(qq1-BS[1][0])+( qq2-BS[1][1])*( qq2-BS[1][1])+(qq5-BS[1][2])*(qq5-BS[1][2]) )-sqrt((qq1-BS[0][0])*(qq1-BS[0][0])+( qq2-BS[0][1])*( qq2-BS[0][1])+(qq5-BS[0][2])*(qq5-BS[0][2] ))-h[0]),2)
  860. +pow(( sqrt((qq1-BS[2][0])*(qq1-BS[2][0])+( qq2-BS[2][1])*( qq2-BS[2][1])+(qq5-BS[2][2])*(qq5-BS[2][2]) )-sqrt((qq1-BS[0][0])*(qq1-BS[0][0])+( qq2-BS[0][1])*( qq2-BS[0][1])+ (qq5-BS[0][2])*(qq5-BS[0][2]))-h[1]),2)
  861. +pow(( sqrt((qq1-BS[3][0])*(qq1-BS[3][0])+( qq2-BS[3][1])*( qq2-BS[3][1])+(qq5-BS[3][2])*(qq5-BS[3][2]) )-sqrt((qq1-BS[0][0])*(qq1-BS[0][0])+( qq2-BS[0][1])*( qq2-BS[0][1])+ (qq5-BS[0][2])*(qq5-BS[0][2]))-h[2]),2);
  862. fitness[0][i]=re;
  863. }
  864. /*产生初始粒子和速度*/
  865. /*个体极值和群体极值*/
  866. double t;double fitnesscopy[1][100] = {0};
  867. for(i=0;i<sizepop;i++)
  868. {
  869. fitnesscopy[0][i] = fitness[0][i];
  870. }
  871. for(i=0;i<sizepop;i++)
  872. {
  873. for(j=i+1;j<sizepop;j++)
  874. {
  875. if (fitnesscopy[0][i]>fitnesscopy[0][j])
  876. {
  877. t = fitnesscopy[0][i];
  878. fitnesscopy[0][i]=fitnesscopy[0][j];
  879. fitnesscopy[0][j] = t;
  880. }
  881. }
  882. }//正确
  883. double bestfitness = fitnesscopy[0][0];
  884. int bestindex,kk;
  885. for(i=0;i<sizepop;i++)
  886. {
  887. if(fitness[0][i]==bestfitness)
  888. {
  889. kk = i;
  890. }
  891. }
  892. bestindex = kk;//正确
  893. double zbest[1][3] = {0};
  894. for (j = 0;j<3;j++)
  895. {
  896. zbest[0][j] = pop[bestindex][j];//全局最佳
  897. }
  898. double gbest[100][3] = {0};
  899. for(i = 0;i<100;i++)
  900. {
  901. for (j = 0;j<3;j++)
  902. {
  903. gbest[i][j] = pop[i][j];//个体最佳
  904. }
  905. }
  906. double fitnessgbest[1][100] = {0};
  907. double fitnesszbest;
  908. for(i = 0;i<100;i++)
  909. {
  910. fitnessgbest[0][i]=fitness[0][i];
  911. }
  912. fitnesszbest = bestfitness;
  913. /*个体极值和群体极值*/
  914. /*迭代寻优*/
  915. double w;
  916. double ws = 0.9;
  917. double we = 0.4;
  918. //srand((unsigned)time(NULL));
  919. for(i = 0;i<maxgen;i++)
  920. {
  921. w = ws - (ws-we)*(i/maxgen)*(i/maxgen);
  922. for(j=0;j<sizepop;j++)
  923. {
  924. //V[j][0]= w*V[j][0]+c1*(rand()+1)/(RAND_MAX+2.0)*(gbest[j][0]-pop[j][0])+c2*(rand()+1)/(RAND_MAX+2.0)*(zbest[j][0]-pop[j][0]);
  925. //V[j][1]= w*V[j][1]+c1*(rand()+1)/(RAND_MAX+2.0)*(gbest[j][1]-pop[j][1])+c2*(rand()+1)/(RAND_MAX+2.0)*(zbest[j][1]-pop[j][1]);
  926. //V[j][2]= w*V[j][2]+c1*(rand()+1)/(RAND_MAX+2.0)*(gbest[j][2]-pop[j][2])+c2*(rand()+1)/(RAND_MAX+2.0)*(zbest[j][2]-pop[j][2]);
  927. V[j][0]= w*V[j][0]+c1*((rand())/double(RAND_MAX))*(gbest[j][0]-pop[j][0])+c2*((rand())/double(RAND_MAX))*(zbest[j][0]-pop[j][0]);//rand(0-1)
  928. V[j][1]= w*V[j][1]+c1*((rand())/double(RAND_MAX))*(gbest[j][1]-pop[j][1])+c2*((rand())/double(RAND_MAX))*(zbest[j][1]-pop[j][1]);
  929. V[j][2]= w*V[j][2]+c1*((rand())/double(RAND_MAX))*(gbest[j][2]-pop[j][2])+c2*((rand())/double(RAND_MAX))*(zbest[j][2]-pop[j][2]);
  930. if(V[j][0]>Vmax)
  931. {
  932. V[j][0]=Vmax;
  933. }
  934. if(V[j][1]>Vmax)
  935. {
  936. V[j][1]=Vmax;
  937. }
  938. if(V[j][2]>Vmax)
  939. {
  940. V[j][2]=Vmax;
  941. }
  942. if(V[j][0]<Vmin)
  943. {
  944. V[j][0]=Vmin;
  945. }
  946. if(V[j][1]<Vmin)
  947. {
  948. V[j][1]=Vmin;
  949. }
  950. if(V[j][2]<Vmin)
  951. {
  952. V[j][2]=Vmin;
  953. }
  954. //种群更新//
  955. pop[j][0]=pop[j][0]+ V[j][0];
  956. pop[j][1]= pop[j][1]+V[j][1];
  957. pop[j][2]= pop[j][2]+V[j][2];
  958. if(pop[j][0]>popmax)
  959. {
  960. pop[j][0]=popmax;
  961. }
  962. if(pop[j][1]>popmax)
  963. {
  964. pop[j][1]=popmax;
  965. }
  966. if(pop[j][2]>popmax)
  967. {
  968. pop[j][2]=popmax;
  969. }
  970. if(pop[j][0]<popmin)
  971. {
  972. pop[j][0]=popmin;
  973. }
  974. if(pop[j][1]<popmin)
  975. {
  976. pop[j][1]=popmin;
  977. }
  978. if(pop[j][2]<popmin)
  979. {
  980. pop[j][2]=popmin;
  981. }
  982. //种群更新//
  983. double qq3,qq4,qq6,re1;
  984. qq3 = pop[j][0];
  985. qq4 = pop[j][1];
  986. qq6= pop[j][2];
  987. re1=pow((sqrt((qq3-BS[1][0])*(qq3-BS[1][0])+( qq4-BS[1][1])*( qq4-BS[1][1])+(qq6-BS[1][2])*(qq6-BS[1][2]) )-sqrt((qq3-BS[0][0])*(qq3-BS[0][0])+( qq4-BS[0][1])*( qq4-BS[0][1])+(qq6-BS[0][2])*(qq6-BS[0][2] ))-h[0]),2)
  988. +pow(( sqrt((qq3-BS[2][0])*(qq3-BS[2][0])+( qq4-BS[2][1])*( qq4-BS[2][1])+(qq6-BS[2][2])*(qq6-BS[2][2]) )-sqrt((qq3-BS[0][0])*(qq3-BS[0][0])+( qq4-BS[0][1])*( qq4-BS[0][1])+ (qq6-BS[0][2])*(qq6-BS[0][2]))-h[1]),2)
  989. +pow(( sqrt((qq3-BS[3][0])*(qq3-BS[3][0])+( qq4-BS[3][1])*( qq4-BS[3][1])+(qq6-BS[3][2])*(qq6-BS[3][2]) )-sqrt((qq3-BS[0][0])*(qq3-BS[0][0])+( qq4-BS[0][1])*( qq4-BS[0][1])+ (qq6-BS[0][2])*(qq6-BS[0][2]))-h[2]),2);
  990. fitness[0][j]=re1;
  991. }//正确
  992. for(j=0;j<sizepop;j++)
  993. {
  994. if( fitness[0][j] < fitnessgbest[0][j])
  995. {
  996. gbest[j][0] = pop[j][0];
  997. gbest[j][1] = pop[j][1];
  998. gbest[j][2] = pop[j][2];
  999. fitnessgbest[0][j] = fitness[0][j];
  1000. }//个体最优更新
  1001. if(fitness[0][j] < fitnesszbest)
  1002. {
  1003. zbest[0][0] = pop[j][0];
  1004. zbest[0][1] = pop[j][1];
  1005. zbest[0][2] = pop[j][2];
  1006. fitnesszbest = fitness[0][j];
  1007. }//群体最优更新
  1008. }
  1009. }
  1010. /*迭代寻优*/
  1011. p->posx = zbest[0][0];
  1012. p->posy = zbest[0][1];
  1013. p->posz = zbest[0][2];
  1014. return p;
  1015. }
  1016. /*
  1017. * 通过fang算法计算tdoa 2维定位结果
  1018. *
  1019. * param
  1020. * pRdm 分站数据集
  1021. *
  1022. * return
  1023. * 返回定位结果
  1024. *
  1025. */
  1026. std::shared_ptr<POS> LocateAlgorithm::tdoa_locate_2d_by_fang(std::shared_ptr<ReceiveDataMap> pRdm)
  1027. {
  1028. if (pRdm->size() < 3)
  1029. {
  1030. return nullptr;
  1031. }
  1032. std::shared_ptr<POS> p = std::make_shared<POS>();
  1033. //2D定位
  1034. std::vector<double> vtk; //保存ki
  1035. std::vector<double> vtd; //保存di
  1036. std::vector<_coordinate> vtc; //保存xi,1 yi,1
  1037. std::vector<_point> vtp; //保存三角形的顶点坐标
  1038. std::vector<int> vts;
  1039. vtk.resize(0);
  1040. vtd.resize(0);
  1041. vtc.resize(0);
  1042. vts.resize(0);
  1043. _coordinate tc;
  1044. int i = 0;
  1045. unsigned long long time_stamp = 0; //保存第一次的插值时间戳
  1046. for (ReceiveDataMap::iterator it = pRdm->begin();it != pRdm->end()&&i<3;++it)
  1047. {
  1048. double k = 0 ;
  1049. k = pow(it->second->x,2) + pow(it->second->y,2) ;//+ pow(it->second->z,2)
  1050. vtk.push_back(k);
  1051. if (i == 0 )
  1052. {
  1053. time_stamp = it->second->rec_time_stamp;
  1054. tc.x = it->second->x;
  1055. tc.y = it->second->y;
  1056. //tc.z = it->second->z;
  1057. }
  1058. long long diff_time = it->second->rec_time_stamp - time_stamp;
  1059. double d = 0;
  1060. d = CFunctions::getDistance(diff_time,CFunctions::TDOA);
  1061. if (i>0)
  1062. {
  1063. if (d>0)
  1064. {
  1065. vts.push_back(1);
  1066. }
  1067. else
  1068. {
  1069. vts.push_back(-1);
  1070. }
  1071. }
  1072. vtd.push_back(d);
  1073. _coordinate dtc;
  1074. dtc.x = it->second->x - tc.x;
  1075. dtc.y = it->second->y - tc.y;
  1076. //dtc.z = it->second->z - tc.z;
  1077. vtc.push_back(dtc);
  1078. _point p;
  1079. p.x = it->second->x;
  1080. p.y = it->second->y;
  1081. vtp.push_back(p);
  1082. i++;
  1083. }
  1084. double a[4] = {0};
  1085. double dt = 0;
  1086. dt = vtd[1]*vtc[2].y - vtd[2]*vtc[1].y;
  1087. if (fabs(dt) < 1E-5)
  1088. {
  1089. return nullptr;
  1090. }
  1091. a[0] = (vtd[2]*vtc[1].x - vtd[1]*vtc[2].x)/dt;
  1092. a[1] = (-1)*((vtk[1] - vtk[0] - pow(vtd[1],2))*vtd[2] - (vtk[2] - vtk[0] - pow(vtd[2],2))*vtd[1])*0.5/dt;
  1093. if (fabs(vtd[1]) < 1E-5)
  1094. {
  1095. return nullptr;
  1096. }
  1097. a[2] = 0.5*(vtk[1] - vtk[0] - pow(vtd[1],2) - 2*vtc[1].y*a[1])/vtd[1];
  1098. a[3] = -1*(vtc[1].x + vtc[1].y*a[0])/vtd[1];
  1099. double A = 0, B = 0, C = 0;
  1100. A = pow(a[3],2) - 1 - pow(a[0],2);
  1101. B = 2*(a[2]*a[3] + tc.x + a[0]*(tc.y - a[1]));
  1102. C = pow(a[2],2) - pow(tc.x,2) - pow(tc.y-a[1],2);
  1103. _point pos[MAX_READER_TDOA_PATH_NUMS];
  1104. double delta = 0.0;
  1105. int count = 0;
  1106. delta = pow(B,2) - 4*A*C;
  1107. if (delta > 0)
  1108. {
  1109. pos[0].x = ((-1)*B + sqrt(delta))/(2*A);
  1110. pos[0].y = a[0]*pos[0].x + a[1];
  1111. //TRACE(_T("x1: %.2f , y1: %.2f \r\n"),pos[0].x,pos[0].y);
  1112. pos[1].x = ((-1)*B - sqrt(delta))/(2*A);
  1113. pos[1].y = a[0]*pos[1].x + a[1];
  1114. //TRACE(_T("x2: %.2f , y2: %.2f \r\n"),pos[1].x,pos[1].y);
  1115. count = 2;
  1116. }else{
  1117. return nullptr;
  1118. }
  1119. int idx = -1;
  1120. int points = 0;
  1121. for (int i=0;i<count;i++)
  1122. {
  1123. //if (IsInTriangle(vtp,pos[i]))
  1124. {
  1125. //TRACE(_T("point is in triangle. the index is : %d \r\n"),i);
  1126. bool condition[2] = {false,false};
  1127. double d1 = 0,d2 = 0;
  1128. //对两解进行判断,需要同时满足两个条件
  1129. //1.解到点1和点2的距离差的方向性和之前的参数相同
  1130. //2.解到点1和点3的距离差的方向性和之前的参数相同
  1131. d1 = sqrt(pow(vtp[1].x - pos[i].x,2) + pow(vtp[1].y - pos[i].y,2)) - sqrt(pow(vtp[0].x - pos[i].x,2) + pow(vtp[0].y - pos[i].y,2));
  1132. d2 = sqrt(pow(vtp[2].x - pos[i].x,2) + pow(vtp[2].y - pos[i].y,2)) - sqrt(pow(vtp[0].x - pos[i].x,2) + pow(vtp[0].y - pos[i].y,2));
  1133. if ((d1 < 0&&vts[0]==-1)||(d1>0&&vts[0] == 1))
  1134. {
  1135. condition[0] = true;
  1136. }
  1137. if ((d2 < 0&&vts[1]==-1)||(d2>0&&vts[1] == 1))
  1138. {
  1139. condition[1] =true;
  1140. }
  1141. if (condition[0]&&condition[1])
  1142. {
  1143. idx = i;
  1144. }
  1145. }//else{
  1146. //TRACE(_T("point is not in triangle. the index is : %d \r\n"),i);
  1147. //}
  1148. }
  1149. if (points == 2)
  1150. {
  1151. idx = -1;
  1152. //TRACE(_T("There hava 2 point.\r\n"));
  1153. }
  1154. if (idx != -1)
  1155. {
  1156. //TRACE(_T("the idx is : %d \r\n"),idx);
  1157. p->posx = pos[idx].x;
  1158. p->posy = pos[idx].y;
  1159. }else{
  1160. return nullptr;
  1161. }
  1162. return p;
  1163. }
  1164. std::vector<std::shared_ptr<POS>> LocateAlgorithm::tdoa_locate_2d_by_fang(std::shared_ptr<ReceiveDataUnorderedMap> pRdm)
  1165. {
  1166. std::vector<std::shared_ptr<POS>> vt_pos;
  1167. vt_pos.resize(0);
  1168. if (pRdm->size() < 3)
  1169. {
  1170. return std::move(vt_pos);
  1171. }
  1172. std::shared_ptr<POS> p = std::make_shared<POS>();
  1173. //2D定位
  1174. std::vector<double> vtk; //保存ki
  1175. std::vector<double> vtd; //保存di
  1176. std::vector<_coordinate> vtc; //保存xi,1 yi,1
  1177. std::vector<_point> vtp; //保存三角形的顶点坐标
  1178. std::vector<int> vts;
  1179. vtk.resize(0);
  1180. vtd.resize(0);
  1181. vtc.resize(0);
  1182. vts.resize(0);
  1183. _coordinate tc;
  1184. int i = 0;
  1185. unsigned long long time_stamp = 0; //保存第一次的插值时间戳
  1186. for (auto it = pRdm->begin();it != pRdm->end()&&i<3;++it)
  1187. {
  1188. double k = 0 ;
  1189. k = pow(it->second->x,2) + pow(it->second->y,2) ;//+ pow(it->second->z,2)
  1190. vtk.push_back(k);
  1191. if (i == 0 )
  1192. {
  1193. time_stamp = it->second->rec_time_stamp;
  1194. tc.x = it->second->x;
  1195. tc.y = it->second->y;
  1196. //tc.z = it->second->z;
  1197. }
  1198. long long diff_time = it->second->rec_time_stamp - time_stamp;
  1199. double d = 0;
  1200. d = CFunctions::getDistance(diff_time,CFunctions::TDOA);
  1201. if (i>0)
  1202. {
  1203. if (d>0)
  1204. {
  1205. vts.push_back(1);
  1206. }
  1207. else
  1208. {
  1209. vts.push_back(-1);
  1210. }
  1211. }
  1212. vtd.push_back(d);
  1213. _coordinate dtc;
  1214. dtc.x = it->second->x - tc.x;
  1215. dtc.y = it->second->y - tc.y;
  1216. //dtc.z = it->second->z - tc.z;
  1217. vtc.push_back(dtc);
  1218. _point p;
  1219. p.x = it->second->x;
  1220. p.y = it->second->y;
  1221. vtp.push_back(p);
  1222. i++;
  1223. }
  1224. double a[4] = {0};
  1225. double dt = 0;
  1226. dt = vtd[1]*vtc[2].y - vtd[2]*vtc[1].y;
  1227. if (fabs(dt) < ZERO_PRECISION)
  1228. {
  1229. //return nullptr;
  1230. return std::move(vt_pos);
  1231. }
  1232. a[0] = (vtd[2]*vtc[1].x - vtd[1]*vtc[2].x)/dt;
  1233. a[1] = (-1)*((vtk[1] - vtk[0] - pow(vtd[1],2))*vtd[2] - (vtk[2] - vtk[0] - pow(vtd[2],2))*vtd[1])*0.5/dt;
  1234. if (fabs(vtd[1]) < ZERO_PRECISION)
  1235. {
  1236. //return nullptr;
  1237. return std::move(vt_pos);
  1238. }
  1239. a[2] = 0.5*(vtk[1] - vtk[0] - pow(vtd[1],2) - 2*vtc[1].y*a[1])/vtd[1];
  1240. a[3] = -1*(vtc[1].x + vtc[1].y*a[0])/vtd[1];
  1241. double A = 0, B = 0, C = 0;
  1242. A = pow(a[3],2) - 1 - pow(a[0],2);
  1243. B = 2*(a[2]*a[3] + tc.x + a[0]*(tc.y - a[1]));
  1244. C = pow(a[2],2) - pow(tc.x,2) - pow(tc.y-a[1],2);
  1245. _point pos[MAX_READER_TDOA_PATH_NUMS];
  1246. double delta = 0.0;
  1247. int count = 0;
  1248. delta = pow(B,2) - 4*A*C;
  1249. if (delta > 0)
  1250. {
  1251. pos[0].x = ((-1)*B + sqrt(delta))/(2*A);
  1252. pos[0].y = a[0]*pos[0].x + a[1];
  1253. //TRACE(_T("x1: %.2f , y1: %.2f \r\n"),pos[0].x,pos[0].y);
  1254. pos[1].x = ((-1)*B - sqrt(delta))/(2*A);
  1255. pos[1].y = a[0]*pos[1].x + a[1];
  1256. //TRACE(_T("x2: %.2f , y2: %.2f \r\n"),pos[1].x,pos[1].y);
  1257. count = 2;
  1258. }else{
  1259. //return nullptr;
  1260. return std::move(vt_pos);
  1261. }
  1262. for (auto i = 0;i < count;i++)
  1263. {
  1264. bool condition[2] = {false,false};
  1265. double d1 = 0,d2 = 0;
  1266. d1 = sqrt(pow(vtp[1].x - pos[i].x,2) + pow(vtp[1].y - pos[i].y,2)) - sqrt(pow(vtp[0].x - pos[i].x,2) + pow(vtp[0].y - pos[i].y,2));
  1267. d2 = sqrt(pow(vtp[2].x - pos[i].x,2) + pow(vtp[2].y - pos[i].y,2)) - sqrt(pow(vtp[0].x - pos[i].x,2) + pow(vtp[0].y - pos[i].y,2));
  1268. if ((d1 < 0&&vts[0]==-1)||(d1>0&&vts[0] == 1))
  1269. {
  1270. condition[0] = true;
  1271. }
  1272. if ((d2 < 0&&vts[1]==-1)||(d2>0&&vts[1] == 1))
  1273. {
  1274. condition[1] =true;
  1275. }
  1276. if (condition[0]&&condition[1])
  1277. {
  1278. std::shared_ptr<POS> p = std::make_shared<POS>();
  1279. p->posx = pos[i].x;
  1280. p->posy = pos[i].y;
  1281. vt_pos.push_back(std::move(p));
  1282. }
  1283. }
  1284. for (std::vector<std::shared_ptr<POS>>::iterator it_pos = vt_pos.begin();it_pos != vt_pos.end();)
  1285. {
  1286. _point p;
  1287. p.x = (*it_pos)->posx;
  1288. p.y = (*it_pos)->posy;
  1289. p.z = 0;
  1290. if (!IsInTriangle(vtp,p))
  1291. {
  1292. it_pos = vt_pos.erase(it_pos);
  1293. }else{
  1294. it_pos++;
  1295. }
  1296. }
  1297. ////根据解是否在三角形内部来取舍点
  1298. //int idx = -1;
  1299. //int points = 0;
  1300. //for (int i = 0;i<count;i++)
  1301. //{
  1302. // if (IsInTriangle(vtp,pos[i]))
  1303. // {
  1304. // TRACE(_T("point is in triangle. the index is : %d \r\n"),i);
  1305. // bool condition[2] = {false,false};
  1306. // double d1 = 0,d2 = 0;
  1307. // //对两解进行判断,需要同时满足两个条件
  1308. // //1.解到点1和点2的距离差的方向性和之前的参数相同
  1309. // //2.解到点1和点3的距离差的方向性和之前的参数相同
  1310. // d1 = sqrt(pow(vtp[1].x - pos[i].x,2) + pow(vtp[1].y - pos[i].y,2)) - sqrt(pow(vtp[0].x - pos[i].x,2) + pow(vtp[0].y - pos[i].y,2));
  1311. // d2 = sqrt(pow(vtp[2].x - pos[i].x,2) + pow(vtp[2].y - pos[i].y,2)) - sqrt(pow(vtp[0].x - pos[i].x,2) + pow(vtp[0].y - pos[i].y,2));
  1312. // if ((d1 < 0&&vts[0]==-1)||(d1>0&&vts[0] == 1))
  1313. // {
  1314. // condition[0] = true;
  1315. // }
  1316. // if ((d2 < 0&&vts[1]==-1)||(d2>0&&vts[1] == 1))
  1317. // {
  1318. // condition[1] =true;
  1319. // }
  1320. // if (condition[0]&&condition[1])
  1321. // {
  1322. // idx = i;
  1323. // }
  1324. // }else{
  1325. // TRACE(_T("point is not in triangle. the index is : %d \r\n"),i);
  1326. // }
  1327. //}
  1328. //if (points == 2)
  1329. //{
  1330. // idx = -1;
  1331. // //TRACE(_T("There hava 2 point.\r\n"));
  1332. //}
  1333. //if (idx != -1)
  1334. //{
  1335. // //TRACE(_T("the idx is : %d \r\n"),idx);
  1336. // p->posx = pos[idx].x;
  1337. // p->posy = pos[idx].y;
  1338. //}else{
  1339. // //return nullptr;
  1340. // return std::move(vt_pos);
  1341. //}
  1342. return std::move(vt_pos);;
  1343. }
  1344. /*
  1345. * 计算tdoa中,由于高度误差存在影响,计算去掉高度误差后的新距离差
  1346. *
  1347. * param
  1348. * md 实际测量距离差
  1349. * d 距离差两分站之间的距离
  1350. * h 高度误差
  1351. *
  1352. * return
  1353. * 返回实际的距离差
  1354. *
  1355. */
  1356. double LocateAlgorithm::cal_real_distance_diff(const double& md,const double& d,const double& h,const int& s)
  1357. {
  1358. double A = 0.0,B = 0.0,C = 0.0;
  1359. A = 2*md;
  1360. B = 2*d;
  1361. C = 2*md*pow(h,2) - pow(d,2) + pow(md,2);
  1362. double delta = 0.0;
  1363. delta = pow(B,2) - 4*A*C;
  1364. if (delta > 0)
  1365. {
  1366. /* int sign = 0;
  1367. double d1 = 0.0, d2 = 0.0;
  1368. for (unsigned int i = 1;i < 3;i++)
  1369. {
  1370. d1 = 0.0, d2 = 0.0;
  1371. d1 = (-1*B + pow(-1,i)*sqrt(delta))/(2*A);
  1372. d2 = d - d1;
  1373. if (d1 - d2 > 0)
  1374. {
  1375. sign = 1;
  1376. }else{
  1377. sign = -1;
  1378. }
  1379. if (sign == s)
  1380. {
  1381. break;
  1382. }
  1383. }*/
  1384. }else if (fabs(delta) < ZERO_PRECISION)
  1385. {
  1386. }else{
  1387. return INVALID_COORDINATE;
  1388. }
  1389. return INVALID_COORDINATE;
  1390. }
  1391. /*
  1392. * TDOA算法实现
  1393. * 函数名:LocatePos,此函数与Pos函数的区别是本函数两两遍历求坐标
  1394. * 遍历求解
  1395. * param
  1396. * pRdm ------ 存放参与计算的分站信息,主要信息包含同步后的时间戳
  1397. * trpm ------ 地图集
  1398. *
  1399. * return
  1400. * 返回最终结果坐标
  1401. *
  1402. */
  1403. std::unique_ptr<POS> LocateAlgorithm::LocatePos(std::shared_ptr<ReceiveDataMap> pRdm,std::shared_ptr<TDOAReaderPathMap> trpm)
  1404. {
  1405. std::unique_ptr<POS> pos(new POS);
  1406. int nNoReaderPathIdx = 0;
  1407. int nFirstReader[MAX_READER_TDOA_PATH_NUMS] = {-1};
  1408. int nSecondReader[MAX_READER_TDOA_PATH_NUMS] ={-1};
  1409. _coordinate r;
  1410. int r_idx = 0;
  1411. r.reason = 0;
  1412. for (ReceiveDataMap::iterator first = pRdm->begin();first != pRdm->end();++first)
  1413. {
  1414. //如果两级都能找到才运行继续后续操作,否则,表明没有此路径地图集
  1415. TDOAReaderPathMap::iterator rdm_it = trpm->find(first->second->reader_id);
  1416. if(rdm_it == trpm->end()){
  1417. continue;
  1418. }
  1419. double ref_dist = 0.0; //两分站之间无路径的距离差
  1420. _coordinate res[MAX_READER_TDOA_PATH_NUMS];
  1421. int res_idx = 0;
  1422. for(int i = 0;i < MAX_READER_TDOA_PATH_NUMS;i++){
  1423. res[i].x = INVALID_COORDINATE;
  1424. res[i].y = INVALID_COORDINATE;
  1425. res[i].z = INVALID_COORDINATE;
  1426. }
  1427. //存储无路径的两分站的id和坐标
  1428. ReceiveData tmp_reader[MAX_READER_TDOA_PATH_NUMS];
  1429. //存储和第一条分站存在路径的分站信息
  1430. ReceiveData tmp_dist_reader[MAX_READER_TDOA_PATH_NUMS];
  1431. //和第一个分站存在地图集的分站个数,在丁字或十字路口,数量可能为3个或4个
  1432. int nDistReaders = 0;
  1433. //获取第一个时间戳
  1434. ReceiveData f1;
  1435. f1.antenna_id = first->second->antenna_id;
  1436. f1.reader_id = first->second->reader_id;
  1437. f1.rec_time_stamp = first->second->rec_time_stamp;
  1438. f1.x = first->second->x;
  1439. f1.y = first->second->y;
  1440. f1.z = first->second->z;
  1441. ReceiveDataMap::iterator second = first;
  1442. std::advance(second,1);
  1443. //从第二个开始遍历
  1444. for(;second != pRdm->end();++second){
  1445. //获取第二个时间戳
  1446. ReceiveData f2;
  1447. f2.antenna_id = second->second->antenna_id;
  1448. f2.reader_id = second->second->reader_id;
  1449. f2.rec_time_stamp = second->second->rec_time_stamp;
  1450. f2.x = second->second->x;
  1451. f2.y = second->second->y;
  1452. f2.z = second->second->z;
  1453. //时间戳异常
  1454. if(f1.rec_time_stamp == LLONG_MAX || f2.rec_time_stamp == LLONG_MAX){
  1455. continue;
  1456. }
  1457. if(f1.reader_id == f2.reader_id){
  1458. continue;
  1459. }
  1460. ReaderPathMap::iterator rpm_it = trpm->find(f1.reader_id)->second->find(f2.reader_id);
  1461. if(rpm_it == trpm->find(f1.reader_id)->second->end()){
  1462. continue;
  1463. }
  1464. //根据距离的正负,后续判断计算位置取舍时使用
  1465. int nSign = 1;
  1466. long long diffTime = f1.rec_time_stamp - f2.rec_time_stamp;
  1467. //计算位置
  1468. double dist = CFunctions::getDistance(diffTime,CFunctions::TDOA);
  1469. double readers_dist = sqrt(pow(f1.x - f2.x,2) + pow(f1.y - f2.y,2));
  1470. if(fabs(dist) - readers_dist > 0){
  1471. continue;
  1472. }
  1473. //如果和第一条分站存在地图集
  1474. tmp_dist_reader[nDistReaders].reader_id = f2.reader_id;
  1475. tmp_dist_reader[nDistReaders].x = f2.x;
  1476. tmp_dist_reader[nDistReaders].y = f2.y;
  1477. tmp_dist_reader[nDistReaders].z = f2.z;
  1478. tmp_dist_reader[nDistReaders].rec_time_stamp = f2.rec_time_stamp;
  1479. nDistReaders++;
  1480. std::shared_ptr<ReaderPath> pRP = trpm->find(f1.reader_id)->second->find(f2.reader_id)->second;
  1481. tmp_reader[nNoReaderPathIdx].reader_id = f2.reader_id;
  1482. tmp_reader[nNoReaderPathIdx].x = f2.x;
  1483. tmp_reader[nNoReaderPathIdx].y = f2.y;
  1484. tmp_reader[nNoReaderPathIdx].z = f2.z;
  1485. tmp_reader[nNoReaderPathIdx].rec_time_stamp = f2.rec_time_stamp;
  1486. nNoReaderPathIdx++;
  1487. //两分站之间的线段个数
  1488. int seg_num = pRP->nRealCalcPoints - 1;
  1489. if(seg_num == 0 || seg_num > 100){
  1490. continue;
  1491. }
  1492. //因为双曲线与分站之间第i条线段或者第j条线段分别有两焦点
  1493. //或者分站之间就一条直线,有两焦点
  1494. double xcross[2] = {INVALID_COORDINATE,INVALID_COORDINATE};
  1495. double ycross[2] = {INVALID_COORDINATE,INVALID_COORDINATE};
  1496. double zcross[2] = {INVALID_COORDINATE,INVALID_COORDINATE};
  1497. int nIdx = 0;
  1498. //根据线段个数开始计算
  1499. for(int i = 0;i < seg_num;i++){
  1500. //计算位置坐标,双曲线和线段相交的交点
  1501. std::shared_ptr<SOLUTION> r = LocateAlgorithm::GetPos(pRP,dist,i);
  1502. //无解或解无效
  1503. if(r == nullptr ){
  1504. continue;
  1505. }
  1506. if(r->nCount == 0){
  1507. continue;
  1508. }
  1509. if(r->nCount == 1){
  1510. xcross[nIdx] = r->x[0];
  1511. ycross[nIdx] = r->y[0];
  1512. res[res_idx].x = xcross[nIdx];
  1513. res[res_idx].y = ycross[nIdx];
  1514. res[res_idx].z = zcross[nIdx];
  1515. nIdx++;
  1516. res_idx++;
  1517. }
  1518. if(r->nCount == 2){
  1519. for(int j = 0; j < 2;j++){
  1520. xcross[j] = r->x[j];
  1521. ycross[j] = r->y[j];
  1522. nIdx++;
  1523. }
  1524. }
  1525. if(nIdx == 2){
  1526. //解到两焦点之间的距离
  1527. double deltad[2] = {0};
  1528. for(int j = 0; j < 2;j ++){
  1529. double d[2] = {0};
  1530. double dx1 = xcross[j] - pRP->x[0];
  1531. double dy1 = ycross[j] - pRP->y[0];
  1532. d[0] = sqrt(pow(dx1,2) + pow(dy1,2));
  1533. double dx2 = xcross[j] - pRP->x[1];
  1534. double dy2 = ycross[j] - pRP->y[1];
  1535. d[1] = sqrt(pow(dx2,2) + pow(dy2,2));
  1536. deltad[j] = d[0] - d[1];
  1537. }
  1538. int idx = 0;
  1539. //应该和之前计算的dist同方向
  1540. if(dist > 0){
  1541. for(int j = 0;j < 2;j++){
  1542. if(deltad[j] > 0){
  1543. idx = j;
  1544. break;
  1545. }
  1546. }
  1547. }else{
  1548. for(int j = 0;j < 2;j++){
  1549. if(deltad[j] < 0){
  1550. idx = j;
  1551. break;
  1552. }
  1553. }
  1554. }
  1555. nFirstReader[res_idx] = f1.reader_id;
  1556. nSecondReader[res_idx] = f2.reader_id;
  1557. res[res_idx].x = xcross[idx];
  1558. res[res_idx].y = ycross[idx];
  1559. res[res_idx].z = zcross[idx];
  1560. res_idx++;
  1561. }
  1562. continue;
  1563. }
  1564. }
  1565. if(res_idx == 1){
  1566. bool bValid = true;
  1567. //如果定位坐标在分站附近,需要判断此分站是否为special,如果为special则接受此定位结果
  1568. for(ReceiveDataMap::iterator it = pRdm->begin();it!=pRdm->end();++it){
  1569. double dist = sqrt(pow(res[0].x - it->second->x,2) + pow(res[0].y - it->second->y,2));
  1570. if(dist<4){
  1571. if(it->second->special == 0){
  1572. r.reason = ALGO_FAILED_CONDITION_12;
  1573. ALGORITHM_FAILED(ALGO_FAILED_CONDITION_12);
  1574. bValid = false;
  1575. }
  1576. }
  1577. }
  1578. if(pRdm->size()>=3&&bValid){
  1579. double dist[2] = {0};
  1580. //更改选取分站方式:
  1581. //选取和第一个分站有地图集的两个分站来参与判别,
  1582. //例如:分站顺序是2,1,3,其中(1,2)和(2,3)之间有地图集,那么就选取1,3来参与判别
  1583. //如果分站顺序是1,2,3,其中(1,2)有地图集,(1,3)之间无地图集,那么此时就不使用以下方法进行判别
  1584. if(nDistReaders == 2){
  1585. dist[0] = sqrt(pow(res[0].x - tmp_dist_reader[0].x,2) + pow(res[0].y - tmp_dist_reader[0].y,2));
  1586. dist[1] = sqrt(pow(res[0].x - tmp_dist_reader[1].x,2) + pow(res[0].y - tmp_dist_reader[1].y,2));
  1587. double dif1 = fabs(dist[0] - dist[1]);
  1588. long long dt = tmp_dist_reader[0].rec_time_stamp - tmp_dist_reader[1].rec_time_stamp;
  1589. double dif2 = fabs(CFunctions::getDistance(dt,CFunctions::TDOA));
  1590. double dif = fabs(dif1 - dif2);
  1591. if(dif>10){
  1592. r.reason = ALGO_FAILED_CONDITION_13;
  1593. ALGORITHM_FAILED(ALGO_FAILED_CONDITION_13);
  1594. bValid = false;
  1595. }
  1596. }
  1597. }
  1598. if(bValid){
  1599. r.x = res[0].x;
  1600. r.y = res[0].y;
  1601. r.z = res[0].z;
  1602. for (int i=0;i<MAX_READER_TDOA_PATH_NUMS;i++)
  1603. {
  1604. pos->dDiff[i] = 0;
  1605. }
  1606. }
  1607. }else if(res_idx >= 2){
  1608. double d = 99999.9;
  1609. long long dt = tmp_reader[0].rec_time_stamp - tmp_reader[1].rec_time_stamp;
  1610. //计算位置
  1611. double ref_dist = abs(CFunctions::getDistance(dt,CFunctions::TDOA));
  1612. for(int i = 0;i<res_idx;i++){
  1613. double d_tmp = abs(sqrt(pow(res[i].x - tmp_reader[0].x,2) + pow(res[i].y - tmp_reader[0].y,2)) - sqrt(pow(res[i].x - tmp_reader[1].x,2) + pow(res[i].y - tmp_reader[1].y,2)));
  1614. double d_diff = abs(d_tmp - ref_dist);
  1615. pos->dDiff[i] = d_diff;
  1616. if(d_diff < d){
  1617. d = d_diff;
  1618. r_idx = i;
  1619. }
  1620. }
  1621. if(pRdm->size()>=3){
  1622. if(pos->dDiff[r_idx] < 4){
  1623. r.x = res[r_idx].x;
  1624. r.y = res[r_idx].y;
  1625. r.z = res[r_idx].z;
  1626. }else{
  1627. r.reason = ALGO_FAILED_CONDITION_14;
  1628. ALGORITHM_FAILED(ALGO_FAILED_CONDITION_14);
  1629. }
  1630. }
  1631. }else{
  1632. //解的个数<=0
  1633. continue;
  1634. }
  1635. //表示找到解了,后面的数据都不可信
  1636. if (r.x != INVALID_COORDINATE && r.y != INVALID_COORDINATE)
  1637. {
  1638. break;
  1639. }
  1640. }
  1641. //最终得到的结果解
  1642. pos->posx = r.x;
  1643. pos->posy = r.y;
  1644. pos->posz = 0;
  1645. pos->reason = r.reason;
  1646. if(pos->posx != INVALID_COORDINATE && pos->posy != INVALID_COORDINATE){
  1647. pos->nFirstReader = nFirstReader[r_idx];
  1648. pos->nSecondReader = nSecondReader[r_idx];
  1649. }
  1650. return pos;
  1651. }
  1652. SOLUTION* LocateAlgorithm::GetPos(ReaderPathMap rpm,int sta_num,int ant,double dist,int i)
  1653. {
  1654. double x1 = 0.0;
  1655. double y1 = 0.0;
  1656. double x2 = 0.0;
  1657. double y2 = 0.0;
  1658. double deta = 0.0; //b^2 - 4*a*c;
  1659. int count = 0; //方程解的个数
  1660. if((rpm.find(sta_num)->second->px[i] - rpm.find(sta_num)->second->px[i+1]) == 0){
  1661. //方程一的解
  1662. x1 = rpm.find(sta_num)->second->px[i];
  1663. y1 = rpm.find(sta_num)->second->y[ant] - dist;
  1664. x2 = rpm.find(sta_num)->second->px[i];
  1665. y2 = dist + rpm.find(sta_num)->second->y[ant];
  1666. count = 2;
  1667. }else{
  1668. //计算斜率
  1669. double k = (rpm.find(sta_num)->second->py[i] - rpm.find(sta_num)->second->py[i+1])/(rpm.find(sta_num)->second->px[i] - rpm.find(sta_num)->second->px[i+1]);
  1670. //方程二
  1671. //确保有两解的情况
  1672. //圆心为s[sta_num].x[ant],s[sta_num].y[ant];
  1673. double a = pow(k,2) + 1;
  1674. double b = 2*(k*(rpm.find(sta_num)->second->py[i] - k*rpm.find(sta_num)->second->px[i] - rpm.find(sta_num)->second->y[ant])- rpm.find(sta_num)->second->x[ant]);
  1675. double d = rpm.find(sta_num)->second->py[i] - k * rpm.find(sta_num)->second->px[i] - rpm.find(sta_num)->second->y[ant];
  1676. double c = pow(rpm.find(sta_num)->second->x[ant],2) + pow(d,2) - pow(dist,2);
  1677. //根据deta = b^2 - 4*a*c 判断方程是否有解
  1678. deta = pow(b,2) - 4*a*c; //
  1679. if(deta>0){
  1680. count = 2;
  1681. }else if(deta == 0){
  1682. count = 1;
  1683. }
  1684. else{
  1685. count = 0;
  1686. }
  1687. if(count == 0){
  1688. return NULL;
  1689. }
  1690. //得两解
  1691. x1 = -(b + sqrt(deta))/(2*a);
  1692. x2 = (-b + sqrt(deta) )/(2*a);
  1693. y1 = k*x1 + rpm.find(sta_num)->second->py[i] - k*rpm.find(sta_num)->second->px[i];
  1694. y2 = k*x2 + rpm.find(sta_num)->second->py[i] - k*rpm.find(sta_num)->second->px[i];
  1695. }
  1696. double x[3] = {0,0,0};
  1697. double y[3] = {0,0,0};
  1698. x[0] = x1;
  1699. x[1] = x2;
  1700. y[0] = y1;
  1701. y[1] = y2;
  1702. if(count>0){
  1703. for(int t = 0;t<count;t++){
  1704. if(x[t] < min(rpm.find(sta_num)->second->px[i],rpm.find(sta_num)->second->px[i+1]) - 0.01
  1705. ||x[t] > max(rpm.find(sta_num)->second->px[i],rpm.find(sta_num)->second->px[i+1]) + 0.01
  1706. ||y[t] < min(rpm.find(sta_num)->second->py[i],rpm.find(sta_num)->second->py[i+1])
  1707. ||y[t] > max(rpm.find(sta_num)->second->py[i],rpm.find(sta_num)->second->py[i+1])){
  1708. x[t] = 0;
  1709. y[t] = 0;
  1710. }
  1711. }
  1712. }
  1713. SOLUTION* result = new SOLUTION;
  1714. memset(result->x,0,3*sizeof(double));
  1715. memset(result->y,0,3*sizeof(double));
  1716. memset(result->z,0,3*sizeof(double));
  1717. result->x[0] = x[0];
  1718. result->x[1] = x[1];
  1719. result->y[0] = y[0];
  1720. result->y[1] = y[1];
  1721. result->z[0] = 0;
  1722. result->z[1] = 0;
  1723. return result;
  1724. }
  1725. /*
  1726. *
  1727. * param
  1728. * trpm
  1729. * sta_num
  1730. * ant
  1731. * dist
  1732. * i
  1733. *
  1734. * return
  1735. *
  1736. */
  1737. std::shared_ptr<SOLUTION> LocateAlgorithm::GetPos(std::shared_ptr<TOFReaderPathMap> trpm,int reader_id,int antenna_idx,double dist,int seg_idx)
  1738. {
  1739. std::shared_ptr<SOLUTION> s = std::make_shared<SOLUTION>();
  1740. double x1 = 0.0;
  1741. double y1 = 0.0;
  1742. double x2 = 0.0;
  1743. double y2 = 0.0;
  1744. double deta = 0.0; //b^2 - 4*a*c;
  1745. int count = 0; //方程解的个数
  1746. if(abs(trpm->find(reader_id)->second->px[seg_idx] - trpm->find(reader_id)->second->px[seg_idx + 1]) < ZERO_PRECISION){
  1747. //方程一的解
  1748. x1 = trpm->find(reader_id)->second->px[seg_idx];
  1749. y1 = trpm->find(reader_id)->second->y[antenna_idx] - dist;
  1750. x2 = trpm->find(reader_id)->second->px[seg_idx];
  1751. y2 = dist + trpm->find(reader_id)->second->y[antenna_idx];
  1752. count = 2;
  1753. }/*else if(abs(trpm->find(sta_num)->second->py[i] - trpm->find(sta_num)->second->py[i+1]) < ZERO_PRECISION){
  1754. y1 = y2 = trpm->find(sta_num)->second->py[i];
  1755. x1 = trpm->find(sta_num)->second->y[ant] - dist;
  1756. x2 = trpm->find(sta_num)->second->y[ant] + dist;
  1757. count = 2;
  1758. }*/else{
  1759. //计算斜率
  1760. double k = (trpm->find(reader_id)->second->py[seg_idx] - trpm->find(reader_id)->second->py[seg_idx + 1])/(trpm->find(reader_id)->second->px[seg_idx] - trpm->find(reader_id)->second->px[seg_idx + 1]);
  1761. //double k = (trpm->find(reader_id)->second->y[seg_idx] - trpm->find(reader_id)->second->y[seg_idx + 1])/(trpm->find(reader_id)->second->x[seg_idx] - trpm->find(reader_id)->second->x[seg_idx + 1]);
  1762. //3.4286
  1763. //double b = trpm->find(reader_id)->second->y[antenna_idx] - k*trpm->find(reader_id)->second->x[antenna_idx];
  1764. double b = trpm->find(reader_id)->second->py[seg_idx] - k*trpm->find(reader_id)->second->px[seg_idx];
  1765. double A = pow(k,2) + 1;
  1766. double B = 2*k*(b - trpm->find(reader_id)->second->y[antenna_idx]) - 2*trpm->find(reader_id)->second->x[antenna_idx];
  1767. double C = pow(trpm->find(reader_id)->second->x[antenna_idx],2) + pow(b - trpm->find(reader_id)->second->y[antenna_idx],2) - pow(dist,2);
  1768. deta = pow(B,2) - 4*A*C;
  1769. //方程二
  1770. //确保有两解的情况
  1771. //圆心为s[sta_num].x[ant],s[sta_num].y[ant];
  1772. //double a = pow(k,2) + 1;
  1773. //double b = 2*(k*(trpm->find(sta_num)->second->py[i] - k*trpm->find(sta_num)->second->px[i] - trpm->find(sta_num)->second->y[ant])- trpm->find(sta_num)->second->x[ant]);
  1774. //double d = trpm->find(sta_num)->second->py[i] - k * trpm->find(sta_num)->second->px[i] - trpm->find(sta_num)->second->y[ant];
  1775. ////double b = 2*k*d - 2*trpm->find(sta_num)->second->px[i];
  1776. //double c = pow(trpm->find(sta_num)->second->x[ant],2) + pow(d,2) - pow(dist,2);
  1777. ////根据deta = b^2 - 4*a*c 判断方程是否有解
  1778. //deta = pow(b,2) - 4*a*c; //
  1779. if(deta>0){
  1780. count = 2;
  1781. }else if(deta == 0){
  1782. count = 1;
  1783. }
  1784. else{
  1785. count = 0;
  1786. }
  1787. if(count == 0){
  1788. return NULL;
  1789. }
  1790. //得两解
  1791. x1 = -(B + sqrt(deta))/(2*A);
  1792. x2 = (-B + sqrt(deta) )/(2*A);
  1793. //x1 = -(b + sqrt(deta))/(2*a);
  1794. //x2 = (-b + sqrt(deta) )/(2*a);
  1795. y1 = k*x1 + trpm->find(reader_id)->second->py[seg_idx] - k*trpm->find(reader_id)->second->px[seg_idx];
  1796. y2 = k*x2 + trpm->find(reader_id)->second->py[seg_idx] - k*trpm->find(reader_id)->second->px[seg_idx];
  1797. }
  1798. double x[3] = {0,0,0};
  1799. double y[3] = {0,0,0};
  1800. x[0] = x1;
  1801. x[1] = x2;
  1802. y[0] = y1;
  1803. y[1] = y2;
  1804. int nValidCount = count;
  1805. if(count > 0){
  1806. for(int t = 0;t < count;t++){
  1807. if(x[t] < min(trpm->find(reader_id)->second->px[seg_idx],trpm->find(reader_id)->second->px[seg_idx+1])
  1808. ||x[t] > max(trpm->find(reader_id)->second->px[seg_idx],trpm->find(reader_id)->second->px[seg_idx+1])
  1809. ||y[t] < min(trpm->find(reader_id)->second->py[seg_idx],trpm->find(reader_id)->second->py[seg_idx+1])
  1810. ||y[t] > max(trpm->find(reader_id)->second->py[seg_idx],trpm->find(reader_id)->second->py[seg_idx+1])){
  1811. x[t] = INVALID_COORDINATE;
  1812. y[t] = INVALID_COORDINATE;
  1813. nValidCount--;
  1814. }
  1815. }
  1816. }
  1817. s->nCount = nValidCount;
  1818. switch (s->nCount)
  1819. {
  1820. case 1:
  1821. for (int i = 0;i < count;i++)
  1822. {
  1823. if (x[i] != INVALID_COORDINATE && y[i] != INVALID_COORDINATE)
  1824. {
  1825. s->x[0] = x[i];
  1826. s->y[0] = y[i];
  1827. break;
  1828. }
  1829. }
  1830. break;
  1831. case 2:
  1832. for (int i = 0;i < s->nCount;i++)
  1833. {
  1834. s->x[i] = x[i];
  1835. s->y[i] = y[i];
  1836. }
  1837. break;
  1838. default:
  1839. break;
  1840. }
  1841. return s;
  1842. }
  1843. /*
  1844. * Tof二维定位实现
  1845. * 采用最小二乘法,对多个方程求最优解
  1846. *
  1847. * param
  1848. * pRdtm 定位数据
  1849. *
  1850. * return
  1851. * 定位结果
  1852. */
  1853. std::shared_ptr<algorithm::POS> LocateAlgorithm::tof_locate_2d(std::shared_ptr<ReceiveDataTofVector>& pRdtv)
  1854. {
  1855. std::shared_ptr<algorithm::POS> sp_pos = std::make_shared<algorithm::POS>();
  1856. //方法一:牛顿迭代法
  1857. //Eigen::VectorXd x(2);
  1858. //Eigen::MatrixXd data(3,3);
  1859. //int i = 0;
  1860. //for (auto it = pRdtv->begin();it != pRdtv->end();++it,++i)
  1861. //{
  1862. // data(i,0) = (*it)->x;
  1863. // data(i,1) = (*it)->y;
  1864. // //data(i,2) = (*it)->z;
  1865. // data(i,2) = (*it)->distance;
  1866. //}
  1867. //TofFunctor functor(data);
  1868. //Eigen::NumericalDiff<TofFunctor> numDiff(functor);
  1869. //Eigen::LevenbergMarquardt<Eigen::NumericalDiff<TofFunctor>,double> lm(numDiff);
  1870. //lm.parameters.maxfev = 2000;
  1871. //lm.parameters.xtol = 1.0e-10;
  1872. //// 计算后,如果计算出结果则返回该结果
  1873. //int ret = lm.minimize(x);
  1874. ////x=-5.56 y=1.06 实际位置(-5.55,1.151)
  1875. ////有偏差数据的测试位置:(-5.399,1.288)
  1876. //if(ret > 0)
  1877. //{
  1878. // double a = x[0];
  1879. // double b = x[1];
  1880. // double c = x[2];
  1881. // //double d = x[3];
  1882. // sp_pos->posx = x[0];
  1883. // sp_pos->posy = x[1];
  1884. // //sp_pos->posz = x[2];
  1885. //}
  1886. //方法二:此方法可以进入现场数据调试
  1887. Eigen::MatrixXf A(3,4);
  1888. Eigen::VectorXf B(3,1);
  1889. Eigen::MatrixXf X;
  1890. int i = 0;
  1891. for (auto it = pRdtv->begin();it != pRdtv->end();++it,++i)
  1892. {
  1893. A(i,0) = A(i,1) = 1;
  1894. A(i,2) = -2*(*it)->x;
  1895. A(i,3) = -2*(*it)->y;
  1896. B(i,0) = pow((*it)->distance,2) - pow((*it)->x,2) - pow((*it)->y,2);
  1897. }
  1898. X = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B);
  1899. double a = X(0);
  1900. double b = X(1);
  1901. double c = X(2);
  1902. double d = X(3);
  1903. //x=-5.56 y=1.06 实际位置(-5.55,1.151)
  1904. //有偏差数据的测试位置:(-5.399,1.288)
  1905. sp_pos->posx = X(2);
  1906. sp_pos->posy = X(3);
  1907. return std::move(sp_pos);
  1908. }
  1909. std::vector<op::point> LocateAlgorithm::calc_pos(const op::point& circle_center, const double& radius, const op::point* pLine)
  1910. {
  1911. vector<op::point> ret;
  1912. double r = radius;
  1913. op::point p1, p2;
  1914. if (abs(circle_center.x - pLine[0].x) < 1E-4&&abs(circle_center.x - pLine[1].x) < 1E-4) {
  1915. p1.x = pLine[0].x;
  1916. p1.y = circle_center.y + r;
  1917. p2.x = pLine[0].x;
  1918. p2.y = circle_center.y - r;
  1919. }
  1920. else {
  1921. double k = (pLine[1].y - pLine[0].y) / (pLine[1].x - pLine[0].x);
  1922. double b = pLine[1].y - k * pLine[1].x;
  1923. double A = pow(k, 2) + 1;
  1924. double B = 2 * k*(b - circle_center.y) - 2 * circle_center.x;
  1925. double C = pow(circle_center.x, 2) + pow(b - circle_center.y, 2) - pow(r, 2);
  1926. double dt = pow(B, 2) - 4 * A*C;
  1927. if (dt <= 0) {
  1928. return ret;
  1929. }
  1930. //得两解
  1931. p1.x = -(B + sqrt(dt)) / (2 * A);
  1932. p2.x = (-B + sqrt(dt)) / (2 * A);
  1933. p1.y = k * p1.x + b;
  1934. p2.y = k * p2.x + b;
  1935. }
  1936. ret.push_back(p1);
  1937. ret.push_back(p2);
  1938. return ret;
  1939. }
  1940. /*
  1941. * TDOA算法实现
  1942. * 函数名:Pos
  1943. * 遍历求解
  1944. * param
  1945. * pRdm ------ 存放参与计算的分站信息,主要信息包含同步后的时间戳
  1946. * trpm ------ 地图集
  1947. *
  1948. * return
  1949. * 返回最终结果坐标
  1950. *
  1951. */
  1952. std::shared_ptr<POS> LocateAlgorithm::Pos(std::shared_ptr<ReceiveDataMap> pRdm,std::shared_ptr<TDOAReaderPathMap> trpm)
  1953. {
  1954. std::shared_ptr<POS> pos = std::make_shared<POS>();
  1955. double a = 0;
  1956. double ref_dist = 0.0; //两分站之间无路径的距离差
  1957. _coordinate res[MAX_READER_TDOA_PATH_NUMS];
  1958. int res_idx = 0;
  1959. for(int i = 0;i < MAX_READER_TDOA_PATH_NUMS;i++){
  1960. res[i].x = INVALID_COORDINATE;
  1961. res[i].y = INVALID_COORDINATE;
  1962. res[i].z = INVALID_COORDINATE;
  1963. }
  1964. bool bFirst = true;
  1965. int nNoReaderPathIdx = 0;
  1966. int nFirstReader[MAX_READER_TDOA_PATH_NUMS] = {-1};
  1967. int nSecondReader[MAX_READER_TDOA_PATH_NUMS] ={-1};
  1968. //存储无路径的两分站的id和坐标
  1969. ReceiveData tmp_reader[MAX_READER_TDOA_PATH_NUMS];
  1970. //存储和第一条分站存在路径的分站信息
  1971. ReceiveData tmp_dist_reader[MAX_READER_TDOA_PATH_NUMS];
  1972. int nDistReaders = 0;
  1973. ReceiveDataMap::iterator first = pRdm->begin();
  1974. ReceiveDataMap::iterator second = first;
  1975. //偏移到第二个元素
  1976. std::advance(second,1);
  1977. //获取第一个时间戳
  1978. ReceiveData f1;
  1979. f1.antenna_id = first->second->antenna_id;
  1980. f1.reader_id = first->second->reader_id;
  1981. f1.rec_time_stamp = first->second->rec_time_stamp;
  1982. f1.x = first->second->x;
  1983. f1.y = first->second->y;
  1984. f1.z = first->second->z;
  1985. //从第二个开始遍历
  1986. for(;second != pRdm->end();++second){
  1987. //获取第二个时间戳
  1988. ReceiveData f2;
  1989. f2.antenna_id = second->second->antenna_id;
  1990. f2.reader_id = second->second->reader_id;
  1991. f2.rec_time_stamp = second->second->rec_time_stamp;
  1992. f2.x = second->second->x;
  1993. f2.y = second->second->y;
  1994. f2.z = second->second->z;
  1995. //时间戳异常
  1996. if(f1.rec_time_stamp == LLONG_MAX || f2.rec_time_stamp == LLONG_MAX){
  1997. continue;
  1998. }
  1999. if(f1.reader_id == f2.reader_id){
  2000. continue;
  2001. }
  2002. //如果两级都能找到才运行继续后续操作,否则,表明没有此路径地图集
  2003. TDOAReaderPathMap::iterator rdm_it = trpm->find(f1.reader_id);
  2004. if(rdm_it == trpm->end()){
  2005. continue;
  2006. }
  2007. ReaderPathMap::iterator rpm_it = trpm->find(f1.reader_id)->second->find(f2.reader_id);
  2008. if(rpm_it == trpm->find(f1.reader_id)->second->end()){
  2009. continue;
  2010. }else{
  2011. //如果和第一条分站存在地图集
  2012. tmp_dist_reader[nDistReaders].reader_id = f2.reader_id;
  2013. tmp_dist_reader[nDistReaders].x = f2.x;
  2014. tmp_dist_reader[nDistReaders].y = f2.y;
  2015. tmp_dist_reader[nDistReaders].z = f2.z;
  2016. tmp_dist_reader[nDistReaders].rec_time_stamp = f2.rec_time_stamp;
  2017. nDistReaders++;
  2018. }
  2019. //根据距离的正负,后续判断计算位置取舍时使用
  2020. int nSign = 1;
  2021. long long diffTime = f1.rec_time_stamp - f2.rec_time_stamp;
  2022. //计算位置
  2023. double dist = CFunctions::getDistance(diffTime,CFunctions::TDOA);
  2024. double readers_dist = sqrt(pow(f1.x - f2.x,2) + pow(f1.y - f2.y,2));
  2025. if(fabs(dist) - readers_dist > 0){
  2026. continue;
  2027. }
  2028. std::shared_ptr<ReaderPath> pRP = trpm->find(f1.reader_id)->second->find(f2.reader_id)->second;
  2029. //如果和第一条分站存在地图集
  2030. tmp_reader[nNoReaderPathIdx].reader_id = f2.reader_id;
  2031. tmp_reader[nNoReaderPathIdx].x = f2.x;
  2032. tmp_reader[nNoReaderPathIdx].y = f2.y;
  2033. tmp_reader[nNoReaderPathIdx].z = f2.z;
  2034. tmp_reader[nNoReaderPathIdx].rec_time_stamp = f2.rec_time_stamp;
  2035. nNoReaderPathIdx++;
  2036. //两分站之间的线段个数
  2037. int seg_num = pRP->nRealCalcPoints - 1;
  2038. if(seg_num == 0 || seg_num > 100){
  2039. continue;
  2040. }
  2041. //因为双曲线与分站之间第i条线段或者第j条线段分别有两焦点
  2042. //或者分站之间就一条直线,有两焦点
  2043. double xcross[2] = {INVALID_COORDINATE,INVALID_COORDINATE};
  2044. double ycross[2] = {INVALID_COORDINATE,INVALID_COORDINATE};
  2045. double zcross[2] = {INVALID_COORDINATE,INVALID_COORDINATE};
  2046. int nIdx = 0;
  2047. //根据线段个数开始计算
  2048. for(int i = 0;i < seg_num;i++){
  2049. //计算位置坐标,双曲线和线段相交的交点
  2050. std::shared_ptr<SOLUTION> r = LocateAlgorithm::GetPos(pRP,dist,i);
  2051. //无解或解无效
  2052. if(r == nullptr ||r->nCount == 0){
  2053. continue;
  2054. }
  2055. if(r->nCount == 1){
  2056. xcross[nIdx] = r->x[0];
  2057. ycross[nIdx] = r->y[0];
  2058. res[res_idx].x = xcross[nIdx];
  2059. res[res_idx].y = ycross[nIdx];
  2060. res[res_idx].z = zcross[nIdx];
  2061. nIdx++;
  2062. res_idx++;
  2063. }
  2064. if(r->nCount == 2){
  2065. for(int j = 0; j < 2;j++){
  2066. xcross[j] = r->x[j];
  2067. ycross[j] = r->y[j];
  2068. nIdx++;
  2069. }
  2070. }
  2071. if(nIdx == 2){
  2072. //解到两焦点之间的距离
  2073. double deltad[2] = {0};
  2074. for(int j = 0; j < 2;j ++){
  2075. double d[2] = {0};
  2076. double dx1 = xcross[j] - pRP->x[0];
  2077. double dy1 = ycross[j] - pRP->y[0];
  2078. d[0] = sqrt(pow(dx1,2) + pow(dy1,2));
  2079. double dx2 = xcross[j] - pRP->x[1];
  2080. double dy2 = ycross[j] - pRP->y[1];
  2081. d[1] = sqrt(pow(dx2,2) + pow(dy2,2));
  2082. deltad[j] = d[0] - d[1];
  2083. }
  2084. int idx = 0;
  2085. //应该和之前计算的dist同方向
  2086. if(dist > 0){
  2087. for(int j = 0;j < 2;j++){
  2088. if(deltad[j] > 0){
  2089. idx = j;
  2090. break;
  2091. }
  2092. }
  2093. }else{
  2094. for(int j = 0;j < 2;j++){
  2095. if(deltad[j] < 0){
  2096. idx = j;
  2097. break;
  2098. }
  2099. }
  2100. }
  2101. nFirstReader[res_idx] = f1.reader_id;
  2102. nSecondReader[res_idx] = f2.reader_id;
  2103. res[res_idx].x = xcross[idx];
  2104. res[res_idx].y = ycross[idx];
  2105. res[res_idx].z = zcross[idx];
  2106. res_idx++;
  2107. }
  2108. continue;
  2109. }
  2110. }
  2111. _coordinate r;
  2112. r.x = INVALID_COORDINATE;
  2113. r.y = INVALID_COORDINATE;
  2114. r.z = INVALID_COORDINATE;
  2115. if(res_idx == 1){
  2116. bool bValid = true;
  2117. //如果定位坐标在分站附近,需要判断此分站是否为special,如果为special则接受此定位结果
  2118. for(ReceiveDataMap::iterator it = pRdm->begin();it!=pRdm->end();++it){
  2119. double dist = sqrt(pow(res[0].x - it->second->x,2) + pow(res[0].y - it->second->y,2));
  2120. if(dist<NEAR_READER){
  2121. if(it->second->special == 0){
  2122. r.reason = ALGO_FAILED_CONDITION_12;
  2123. ALGORITHM_FAILED(ALGO_FAILED_CONDITION_12);
  2124. bValid = false;
  2125. }
  2126. }
  2127. }
  2128. if(pRdm->size()>=3&&bValid){
  2129. double dist[2] = {0};
  2130. //更改选取分站方式:
  2131. //选取和第一个分站有地图集的两个分站来参与判别,
  2132. //例如:分站顺序是2,1,3,其中(1,2)和(2,3)之间有地图集,那么就选取1,3来参与判别
  2133. //如果分站顺序是1,2,3,其中(1,2)有地图集,(1,3)之间无地图集,那么此时就不使用以下方法进行判别
  2134. if(nDistReaders == 2){
  2135. //如果第二个分站与这两个分站无地图集则不该使用如下逻辑
  2136. dist[0] = sqrt(pow(res[0].x - tmp_dist_reader[0].x,2) + pow(res[0].y - tmp_dist_reader[0].y,2));
  2137. dist[1] = sqrt(pow(res[0].x - tmp_dist_reader[1].x,2) + pow(res[0].y - tmp_dist_reader[1].y,2));
  2138. double dif1 = fabs(dist[0] - dist[1]);
  2139. long long dt = tmp_dist_reader[0].rec_time_stamp - tmp_dist_reader[1].rec_time_stamp;
  2140. double dif2 = fabs(CFunctions::getDistance(dt,CFunctions::TDOA));
  2141. double dif = fabs(dif1 - dif2);
  2142. if(dif>MIN_DIFFER_DISTANCE){
  2143. r.reason = ALGO_FAILED_CONDITION_13;
  2144. ALGORITHM_FAILED(ALGO_FAILED_CONDITION_13);
  2145. bValid = false;
  2146. }
  2147. }
  2148. }
  2149. if(bValid){
  2150. r.x = res[0].x;
  2151. r.y = res[0].y;
  2152. r.z = res[0].z;
  2153. for (int i=0;i<MAX_READER_TDOA_PATH_NUMS;i++)
  2154. {
  2155. pos->dDiff[i] = 0;
  2156. }
  2157. }
  2158. }
  2159. int r_idx = 0;
  2160. if(res_idx >= 2){
  2161. double d = 99999.9;
  2162. long long dt = tmp_reader[0].rec_time_stamp - tmp_reader[1].rec_time_stamp;
  2163. //计算位置
  2164. double ref_dist = abs(CFunctions::getDistance(dt,CFunctions::TDOA));
  2165. for(int i = 0;i<res_idx;i++){
  2166. double d_tmp = abs(sqrt(pow(res[i].x - tmp_reader[0].x,2) + pow(res[i].y - tmp_reader[0].y,2)) - sqrt(pow(res[i].x - tmp_reader[1].x,2) + pow(res[i].y - tmp_reader[1].y,2)));
  2167. double d_diff = abs(d_tmp - ref_dist);
  2168. pos->dDiff[i] = d_diff;
  2169. if(d_diff < d){
  2170. d = d_diff;
  2171. r_idx = i;
  2172. }
  2173. }
  2174. if(pRdm->size() >= 3){
  2175. if(pos->dDiff[r_idx] < NEAR_READER){
  2176. r.x = res[r_idx].x;
  2177. r.y = res[r_idx].y;
  2178. r.z = res[r_idx].z;
  2179. }else{
  2180. r.reason = ALGO_FAILED_CONDITION_14;
  2181. ALGORITHM_FAILED(ALGO_FAILED_CONDITION_14);
  2182. }
  2183. }
  2184. }
  2185. //最终得到的结果解
  2186. pos->posx = r.x;
  2187. pos->posy = r.y;
  2188. pos->posz = 0;
  2189. pos->reason = r.reason;
  2190. if(pos->posx != INVALID_COORDINATE && pos->posy != INVALID_COORDINATE){
  2191. pos->nFirstReader = nFirstReader[r_idx];
  2192. pos->nSecondReader = nSecondReader[r_idx];
  2193. }
  2194. return pos;
  2195. }
  2196. std::shared_ptr<POS> LocateAlgorithm::Pos(std::shared_ptr<POS> pos,std::shared_ptr<TDOAReaderPathMap> trpm)
  2197. {
  2198. if(trpm->find(pos->nFirstReader) == trpm->end()){
  2199. return nullptr;
  2200. }
  2201. if(trpm->find(pos->nFirstReader)->second->find(pos->nSecondReader)==trpm->find(pos->nFirstReader)->second->end()){
  2202. return nullptr;
  2203. }
  2204. std::shared_ptr<POS> p = std::make_shared<POS>();
  2205. _coordinate res[2];
  2206. int res_idx = 0;
  2207. for(int i = 0;i < 2;i++){
  2208. res[i].x = INVALID_COORDINATE;
  2209. res[i].y = INVALID_COORDINATE;
  2210. res[i].z = INVALID_COORDINATE;
  2211. }
  2212. //获得地图集
  2213. std::shared_ptr<ReaderPath> pRP = trpm->find(pos->nFirstReader)->second->find(pos->nSecondReader)->second;
  2214. //距离差
  2215. double distance = 0.0;
  2216. distance = sqrt(pow(pos->posx - pRP->x[0],2 ) + pow(pos->posy - pRP->y[0],2 )) -
  2217. sqrt(pow(pos->posx - pRP->x[1],2 ) + pow(pos->posy - pRP->y[1],2 ));
  2218. //两分站之间的线段个数
  2219. int seg_num = pRP->nRealCalcPoints - 1;
  2220. if(seg_num == 0){
  2221. return nullptr;
  2222. }
  2223. double xcross[2] = {INVALID_COORDINATE,INVALID_COORDINATE};
  2224. double ycross[2] = {INVALID_COORDINATE,INVALID_COORDINATE};
  2225. double zcross[2] = {INVALID_COORDINATE,INVALID_COORDINATE};
  2226. int nIdx = 0;
  2227. int nResIdx = 0;
  2228. for(int i = 0;i < seg_num;i++){
  2229. //计算位置坐标,双曲线和线段相交的交点
  2230. std::shared_ptr<SOLUTION> r = GetPos(pRP,distance,i);
  2231. //无解或解无效
  2232. if(r == nullptr){
  2233. continue;
  2234. }
  2235. if(r->nCount == 0){
  2236. continue;
  2237. }
  2238. if(r->nCount == 1){
  2239. xcross[nIdx] = r->x[0];
  2240. ycross[nIdx] = r->y[0];
  2241. nIdx++;
  2242. }
  2243. if(r->nCount == 2){
  2244. for(int j = 0; j < 2;j++){
  2245. xcross[j] = r->x[j];
  2246. ycross[j] = r->y[j];
  2247. nIdx++;
  2248. }
  2249. }
  2250. if(nIdx == 2){
  2251. //解到两焦点之间的距离
  2252. double deltad[2] = {0};
  2253. for(int j = 0; j < 2;j ++){
  2254. double d[2] = {0};
  2255. double dx1 = xcross[j] - pRP->x[0];
  2256. double dy1 = ycross[j] - pRP->y[0];
  2257. d[0] = sqrt(pow(dx1,2) + pow(dy1,2));
  2258. double dx2 = xcross[j] - pRP->x[1];
  2259. double dy2 = ycross[j] - pRP->y[1];
  2260. d[1] = sqrt(pow(dx2,2) + pow(dy2,2));
  2261. deltad[j] = d[0] - d[1];
  2262. }
  2263. int idx = 0;
  2264. //应该和之前计算的dist同方向
  2265. if(distance > 0){
  2266. for(int j = 0;j < 2;j++){
  2267. if(deltad[j] > 0){
  2268. idx = j;
  2269. break;
  2270. }
  2271. }
  2272. }else{
  2273. for(int j = 0;j < 2;j++){
  2274. if(deltad[j] < 0){
  2275. idx = j;
  2276. break;
  2277. }
  2278. }
  2279. }
  2280. res[nResIdx].x = xcross[idx];
  2281. res[nResIdx].y = ycross[idx];
  2282. res[nResIdx].z = zcross[idx];
  2283. }
  2284. continue;
  2285. }
  2286. if (res[nResIdx].x == INVALID_COORDINATE && res[nResIdx].y == INVALID_COORDINATE)
  2287. {
  2288. return nullptr;
  2289. }
  2290. p->posx = res[nResIdx].x;
  2291. p->posy = res[nResIdx].y;
  2292. p->posz = res[nResIdx].z;
  2293. return p;
  2294. }
  2295. /*
  2296. * TDOA算法求解实现
  2297. * 函数名:GetPos
  2298. * 实现直线和双曲线求交点
  2299. * param
  2300. * pRP ------ 分站路径
  2301. * dist ------ 距离差
  2302. * i ------ 第i条线段
  2303. *
  2304. * return
  2305. * 满足条件的解,最多两个
  2306. *
  2307. */
  2308. typedef struct
  2309. {
  2310. double x;
  2311. double y;
  2312. double z;
  2313. }PositionStr;
  2314. int calc_card_position(ReceiveData stationFirst, ReceiveData stationSecond, double distanceDiff, PositionStr& result)
  2315. {
  2316. double station_distance = sqrt(pow((stationFirst.x - stationSecond.x), 2) + pow((stationFirst.y - stationSecond.y),2));
  2317. double diff_check = abs(abs(station_distance) - abs(distanceDiff));
  2318. if(abs(distanceDiff) > station_distance)
  2319. {
  2320. //debug_print_syslog(0, "calc_card_position, distanceDiff is longer than stationDistance, read_id_first: %d,read_id_second: %d, diff_check: %f,\
  2321. // first_time: %I64u, second_time: %I64u, station_distance: %f, distanceDiff: %f",
  2322. // stationFirst.reader_id, stationSecond.reader_id,
  2323. // diff_check, stationFirst.rec_time_stamp, stationSecond.rec_time_stamp, station_distance, distanceDiff);
  2324. return -1;
  2325. }
  2326. if(diff_check < 4)
  2327. {
  2328. //debug_print_syslog(0, "calc_card_position, read_id_first: %d,read_id_second: %d, diff_check: %f,\
  2329. // first_time: %I64u, second_time: %I64u, station_distance: %f, distanceDiff: %f",
  2330. // stationFirst.reader_id, stationSecond.reader_id,
  2331. // diff_check, stationFirst.rec_time_stamp, stationSecond.rec_time_stamp, station_distance, distanceDiff);
  2332. return -1;
  2333. }
  2334. double m = station_distance + distanceDiff;
  2335. double n = station_distance - distanceDiff;
  2336. result.x = ((m * stationSecond.x) + (n * stationFirst.x))/(m + n);
  2337. result.y = ((m * stationSecond.y) + (n * stationFirst.y))/(m + n);
  2338. result.z = 0;
  2339. return 0;
  2340. }
  2341. std::unique_ptr<POS> LocateAlgorithm::CalcCardPosition(std::shared_ptr<ReceiveDataMap> pRdm,std::shared_ptr<TDOAReaderPathMap> trpm)
  2342. {
  2343. std::unique_ptr<POS> pos(new POS);
  2344. double a = 0;
  2345. double ref_dist = 0.0; //两分站之间无路径的距离差
  2346. const int Totals = 10;
  2347. _coordinate res;
  2348. int res_idx = 0;
  2349. bool bFirst = true;
  2350. int nNoReaderPathIdx = 0;
  2351. int nFirstReader[Totals] = {-1};
  2352. int nSecondReader[Totals] ={-1};
  2353. pos->posx = INVALID_COORDINATE;
  2354. pos->posy = INVALID_COORDINATE;
  2355. pos->posz = INVALID_COORDINATE;
  2356. pos->reason = ALGO_FAILED_CONDITION_1;
  2357. //存储无路径的两分站的id和坐标
  2358. ReceiveData tmp_reader[Totals];
  2359. //存储和第一条分站存在路径的分站信息
  2360. ReceiveData tmp_dist_reader[Totals];
  2361. int nDistReaders = 0;
  2362. for(ReceiveDataMap::iterator it_tmp = pRdm->begin(); it_tmp != pRdm->end(); it_tmp++)
  2363. {
  2364. //debug_print_syslog(0, "Calculate data, read_id_second: %d, rec_time_stamp: %I64u",
  2365. // it_tmp->second->reader_id, it_tmp->second->rec_time_stamp);
  2366. }
  2367. ReceiveDataMap::iterator first = pRdm->begin();
  2368. ReceiveDataMap::iterator second = first;
  2369. //偏移到第二个元素
  2370. std::advance(second,1);
  2371. //获取第一个时间戳
  2372. ReceiveData f1;
  2373. f1.antenna_id = first->second->antenna_id;
  2374. f1.reader_id = first->second->reader_id;
  2375. f1.rec_time_stamp = first->second->rec_time_stamp;
  2376. f1.x = first->second->x;
  2377. f1.y = first->second->y;
  2378. f1.z = first->second->z;
  2379. //从第二个开始遍历
  2380. for(; second != pRdm->end(); ++second)
  2381. {
  2382. //获取第二个时间戳
  2383. ReceiveData f2;
  2384. f2.antenna_id = second->second->antenna_id;
  2385. f2.reader_id = second->second->reader_id;
  2386. f2.rec_time_stamp = second->second->rec_time_stamp;
  2387. f2.x = second->second->x;
  2388. f2.y = second->second->y;
  2389. f2.z = second->second->z;
  2390. //时间戳异常
  2391. if(f1.rec_time_stamp == LLONG_MAX || f2.rec_time_stamp == LLONG_MAX){
  2392. continue;
  2393. }
  2394. if(f1.reader_id == f2.reader_id){
  2395. continue;
  2396. }
  2397. //如果两级都能找到才运行继续后续操作,否则,表明没有此路径地图集
  2398. TDOAReaderPathMap::iterator rdm_it = trpm->find(f1.reader_id);
  2399. if(rdm_it == trpm->end()){
  2400. continue;
  2401. }
  2402. ReaderPathMap::iterator rpm_it = trpm->find(f1.reader_id)->second->find(f2.reader_id);
  2403. if(rpm_it == trpm->find(f1.reader_id)->second->end()){
  2404. continue;
  2405. }else{
  2406. //如果和第一条分站存在地图集
  2407. tmp_dist_reader[nDistReaders].reader_id = f2.reader_id;
  2408. tmp_dist_reader[nDistReaders].x = f2.x;
  2409. tmp_dist_reader[nDistReaders].y = f2.y;
  2410. tmp_dist_reader[nDistReaders].z = f2.z;
  2411. tmp_dist_reader[nDistReaders].rec_time_stamp = f2.rec_time_stamp;
  2412. nDistReaders++;
  2413. }
  2414. //根据距离的正负,后续判断计算位置取舍时使用
  2415. int nSign = 1;
  2416. long long diffTime = f1.rec_time_stamp - f2.rec_time_stamp;
  2417. //计算位置
  2418. double distDiff = CFunctions::getDistance(diffTime,CFunctions::TDOA);
  2419. PositionStr stationFirst;
  2420. stationFirst.x = f1.x;
  2421. stationFirst.y = f1.y;
  2422. PositionStr stationSecond;
  2423. stationSecond.x = f2.x;
  2424. stationSecond.y = f2.y;
  2425. PositionStr result;
  2426. int rt_val = 0;
  2427. rt_val = calc_card_position(f1, f2, distDiff, result);
  2428. if(0 == rt_val)
  2429. {
  2430. res.x = result.x;
  2431. res.y = result.y;
  2432. res.z = result.z;
  2433. res.reason = ALGO_LOC_SUCCESSED;
  2434. res_idx++;
  2435. pos->posx = res.x;
  2436. pos->posy = res.y;
  2437. pos->posz = res.z;
  2438. pos->reason = res.reason;
  2439. return pos;
  2440. }
  2441. }
  2442. return pos;
  2443. }
  2444. std::shared_ptr<SOLUTION> LocateAlgorithm::GetPos(std::shared_ptr<ReaderPath> pRP,double dist,int i)
  2445. {
  2446. //解的个数
  2447. int count = 0;
  2448. //解的坐标
  2449. double x[2];
  2450. double y[2];
  2451. //double z[2];
  2452. for (int t=0;t<2;t++)
  2453. {
  2454. x[t] = y[t] = INVALID_COORDINATE;
  2455. }
  2456. //双曲线的两个焦点分别是
  2457. double x1 = pRP->x[0];
  2458. double y1 = pRP->y[0];
  2459. double z1 = pRP->z[0];
  2460. double x2 = pRP->x[1];
  2461. double y2 = pRP->y[1];
  2462. double z2 = pRP->z[1];
  2463. if(pRP->px[i] - pRP->px[i+1] == 0){
  2464. //x相等,双曲线在Y轴上
  2465. x[0] = x[1] = pRP->px[i];
  2466. y[0] = (y1 + y2 + dist)/2.0;
  2467. y[1] = (y1 + y2 - dist)/2.0;
  2468. count = 2;
  2469. }
  2470. else{
  2471. //双曲线和直线相交
  2472. //直线常数求解
  2473. double k = (pRP->py[i+1] - pRP->py[i])/(pRP->px[i+1] - pRP->px[i]);
  2474. double kb = pRP->py[i] - k*pRP->px[i];
  2475. //求解的常数中间量
  2476. double d1 = 2*(x2 - x1);//m
  2477. double d2 = 2*(y2 - y1);//n
  2478. double d3 = pow(y1,2) + pow(x1,2)- pow(dist,2) - pow(y2,2) - pow(x2,2);//p
  2479. //方程ax^2 + bx + c = 0
  2480. double a = pow(d1 + d2*k,2) - 4*pow(dist,2)*(pow(k,2) + 1);
  2481. double b = 2*((d3 + d2*kb)*(d1 + d2*k) - 4*pow(dist,2)*(k*(kb - y2) - x2));
  2482. double c = pow(d3 + d2*kb,2) - 4*pow(dist,2)*(pow(x2,2) + pow(kb-y2,2));
  2483. double delta = pow(b,2) - 4*a*c;
  2484. if(delta > 0){
  2485. count = 2;
  2486. }
  2487. else if(delta == 0){
  2488. count = 1;
  2489. }
  2490. else{
  2491. count = 0;
  2492. }
  2493. if(count == 0){
  2494. return nullptr;
  2495. }
  2496. //计算解
  2497. x[0] = -(b + sqrt(delta))/(2*a);
  2498. x[1] = (-b + sqrt(delta) )/(2*a);
  2499. y[0] = k*x[0] + kb;
  2500. y[1] = k*x[1] + kb;
  2501. }
  2502. int nIdx = 0;
  2503. int nValidCount = count;
  2504. //引入误差
  2505. double deviation = 0.0;
  2506. //判断两解是否在线段范围内
  2507. if(count > 0){
  2508. for(int t = 0; t < count ;t++){
  2509. //两分站之间第i条线段的两端点
  2510. if(x[t] < min(pRP->px[i],pRP->px[i+1]) - deviation||
  2511. x[t] > max(pRP->px[i],pRP->px[i+1]) + deviation||
  2512. y[t] < min(pRP->py[i],pRP->py[i+1]) - deviation||
  2513. y[t] > max(pRP->py[i],pRP->py[i+1]) + deviation
  2514. )
  2515. {
  2516. x[t] = INVALID_COORDINATE;
  2517. y[t] = INVALID_COORDINATE;
  2518. nValidCount--;
  2519. nIdx = t;
  2520. }
  2521. }
  2522. }
  2523. std::shared_ptr<SOLUTION> s = std::make_shared<SOLUTION>();
  2524. s->nCount = nValidCount;
  2525. nIdx = -1;
  2526. switch(nValidCount){
  2527. case 1:
  2528. for(int i = 0;i<2;i++){
  2529. if(x[i] != INVALID_COORDINATE && y[i] != INVALID_COORDINATE){
  2530. nIdx = i;
  2531. }
  2532. }
  2533. s->x[0] = x[nIdx];
  2534. s->y[0] = y[nIdx];
  2535. break;
  2536. case 2:
  2537. s->x[0] = x[0];
  2538. s->x[1] = x[1];
  2539. s->y[0] = y[0];
  2540. s->y[1] = y[1];
  2541. break;
  2542. }
  2543. return s;
  2544. }
  2545. /*
  2546. * 判断坐标在地图集上
  2547. *
  2548. * param
  2549. * pos ------ 定位的结果
  2550. * trpm ------ 地图集
  2551. *
  2552. * return
  2553. * false,不在地图集上,true在地图集上
  2554. *
  2555. */
  2556. bool LocateAlgorithm::IsOnMap(std::shared_ptr<POS>& pos,std::shared_ptr<TDOAReaderPathMap> trpm)
  2557. {
  2558. //如果地图集中都找不到此分站开始的地图集,就返回false
  2559. //if(trpm->find(pos->nFirstReader) == trpm->end()){
  2560. // return false;
  2561. //}
  2562. ////如果地图集找不到此两分站的地图集,就返回false
  2563. //if(trpm->find(pos->nFirstReader)->second->find(pos->nSecondReader) == trpm->find(pos->nFirstReader)->second->end())
  2564. //{
  2565. // return false;
  2566. //}
  2567. ////保存两分站的坐标
  2568. //double x[2] = {0};
  2569. //double y[2] = {0};
  2570. //for(int i=0;i<2;i++){
  2571. // x[i] = trpm->find(pos->nFirstReader)->second->find(pos->nSecondReader)->second->x[i];
  2572. // y[i] = trpm->find(pos->nFirstReader)->second->find(pos->nSecondReader)->second->y[i];
  2573. //}
  2574. //if(x[0] == x[1]){
  2575. // //误差1cm
  2576. // if(abs(pos->posx - x[0]) > 1E-2){
  2577. // return false;
  2578. // }
  2579. //}else{
  2580. // double k = (y[1] - y[0])/(x[1] - x[0]);
  2581. // double b = y[1] - k*x[1];
  2582. // double calc_y = k*pos->posx + b;
  2583. // if(abs(calc_y - pos->posy) > 1E-2){
  2584. // return false;
  2585. // }
  2586. //}
  2587. bool bExist = false;
  2588. for (TDOAReaderPathMap::iterator first = trpm->begin();first!=trpm->end();++first)
  2589. {
  2590. for (ReaderPathMap::iterator second = first->second->begin();second != first->second->end();++second)
  2591. {
  2592. _point p,start_p,end_p;
  2593. p.x = pos->posx;
  2594. p.y = pos->posy;
  2595. start_p.x = second->second->px[0];
  2596. start_p.y = second->second->py[0];
  2597. end_p.x = second->second->px[1];
  2598. end_p.y = second->second->py[1];
  2599. bExist = LocateAlgorithm::IsInLine(p,start_p,end_p);
  2600. if (bExist)
  2601. {
  2602. pos->nFirstReader = first->first;
  2603. pos->nSecondReader = second->first;
  2604. return bExist;
  2605. }
  2606. ////保存两分站的坐标
  2607. //double x[2] = {0};
  2608. //double y[2] = {0};
  2609. //for(int i=0;i<2;i++){
  2610. // x[i] = second->second->x[i];
  2611. // y[i] = second->second->y[i];
  2612. //}
  2613. //if(abs(x[0] - x[1]) < 1E-4 && abs(pos->posx - x[0]) < 1E-4){
  2614. // //误差1cm
  2615. // if ((pos->posy > y[0] && pos->posy < y[1])||(pos->posy > y[1] && pos->posy < y[0]))
  2616. // {
  2617. // bExist = true;
  2618. // return bExist;
  2619. // }
  2620. //}else{
  2621. // double k = (y[1] - y[0])/(x[1] - x[0]);
  2622. // double b = y[1] - k*x[1];
  2623. // double calc_y = k*pos->posx + b;
  2624. // if(abs(calc_y - pos->posy) < 1){
  2625. // bExist = true;
  2626. // return bExist;
  2627. // }
  2628. //}
  2629. }
  2630. }
  2631. return bExist;
  2632. }
  2633. bool LocateAlgorithm::IsInTriangle(std::vector<_point> vtp,_point p)
  2634. {
  2635. double sabc = 0,sadb = 0,sbdc = 0,sadc = 0;
  2636. sabc = GetTriangleArea(vtp[0],vtp[1],vtp[2]);
  2637. sadb = GetTriangleArea(vtp[0],p,vtp[1]);
  2638. sbdc = GetTriangleArea(vtp[1],p,vtp[2]);
  2639. sadc = GetTriangleArea(vtp[0],p,vtp[2]);
  2640. double sum = 0.0;
  2641. sum = sadb + sbdc + sadc;
  2642. if ((sabc - sum) > -1E-5 && (sabc - sum) < 1E-5)
  2643. {
  2644. return true;
  2645. }
  2646. else
  2647. {
  2648. return false;
  2649. }
  2650. }
  2651. double LocateAlgorithm::GetTriangleArea(_point p0,_point p1,_point p2)
  2652. {
  2653. _point ab,bc;
  2654. ab.x = p1.x - p0.x;
  2655. ab.y = p1.y - p0.y;
  2656. bc.x = p2.x - p1.x;
  2657. bc.y = p2.y - p1.y;
  2658. return abs(ab.x*bc.y - ab.y*bc.x)/2.0;
  2659. }
  2660. /*
  2661. * 点P是否在以start_p和end_p的线段上
  2662. *
  2663. * param
  2664. * p 识别点
  2665. * start_p 线段端点
  2666. * end_p 线段终点
  2667. *
  2668. * return
  2669. * 在线段上返回true,否则返回false
  2670. *
  2671. */
  2672. bool LocateAlgorithm::IsInLine(_point p,_point start_p,_point end_p)
  2673. {
  2674. double d1 = 0, d2 = 0, d3 = 0;
  2675. d1 = sqrt(pow(p.x - start_p.x,2) + pow(p.y - start_p.y,2));
  2676. d2 = sqrt(pow(p.x - end_p.x,2) + pow(p.y - end_p.y,2));
  2677. d3 = sqrt(pow(start_p.x - end_p.x,2) + pow(start_p.y - end_p.y,2));
  2678. double d4 = 0;
  2679. d4 = abs(d3 - (d1 + d2));
  2680. if (d4 < OFFSET_THRE_IN_LINE)
  2681. {
  2682. //增加点是否在start_p和end_p为顶点的矩形内
  2683. if (abs(start_p.x - end_p.x) < ZERO_PRECISION)
  2684. {
  2685. if (p.y > min(start_p.y,end_p.y) && p.y < max(start_p.y,end_p.y))
  2686. {
  2687. return true;
  2688. }else{
  2689. return false;
  2690. }
  2691. }else if(abs(start_p.y - end_p.y) < ZERO_PRECISION){
  2692. if(p.x > min(start_p.x,end_p.x) && p.x < max(start_p.x,end_p.x)) {
  2693. return true;
  2694. }else{
  2695. return false;
  2696. }
  2697. }else{
  2698. if(p.x > min(start_p.x,end_p.x) && p.x < max(start_p.x,end_p.x)
  2699. && p.y > min(start_p.y,end_p.y) && p.y < max(start_p.y,end_p.y))
  2700. {
  2701. return true;
  2702. }else{
  2703. return false;
  2704. }
  2705. }
  2706. }
  2707. return false;
  2708. }
  2709. bool LocateAlgorithm::PointIsInRect(_point p,_point tp,_point lp,_point bp,_point rp)
  2710. {
  2711. //任意四边形有4个顶点
  2712. int nCount = 4;
  2713. _point RectPoints[4] = { tp, lp, bp, rp };
  2714. int nCross = 0;
  2715. for (int i = 0; i < nCount; i++)
  2716. {
  2717. //依次取相邻的两个点
  2718. _point pStart = RectPoints[i];
  2719. _point pEnd = RectPoints[(i + 1) % nCount];
  2720. //相邻的两个点是平行于x轴的,肯定不相交,忽略
  2721. if ( pStart.y == pEnd.y )
  2722. continue;
  2723. //交点在pStart,pEnd的延长线上,pCur肯定不会与pStart.pEnd相交,忽略
  2724. if ( p.y < min(pStart.y, pEnd.y) || p.y > max( pStart.y, pEnd.y ) )
  2725. continue;
  2726. //求当前点和x轴的平行线与pStart,pEnd直线的交点的x坐标
  2727. double x = (double)(p.y - pStart.y) * (double)(pEnd.x - pStart.x) / (double)(pEnd.y - pStart.y) + pStart.x;
  2728. //若x坐标大于当前点的坐标,则有交点
  2729. if ( x > p.x )
  2730. nCross++;
  2731. }
  2732. // 单边交点为偶数,点在多边形之外
  2733. return (nCross % 2 == 1);
  2734. }
  2735. int LocateAlgorithm::CalcTdoaPosition(std::shared_ptr<ReceiveDataMap> pRdm,std::shared_ptr<TDOAReaderPathMap> trpm,std::vector<std::shared_ptr<POS>>& udm_pos)
  2736. {
  2737. //这部分应该放到组装数据后检测里
  2738. //一:选择基准分站,条件如下:
  2739. //1.基准分站选择:当前分站和后一个分站有地图集,则选为基准分站
  2740. //2.当前分站和后续所有分站无地图集,则当前分站不是基准分站
  2741. //3.重复1,2过程
  2742. for (ReceiveDataMap::iterator first = pRdm->begin();first != pRdm->end();)
  2743. {
  2744. ReceiveDataMap::iterator second = first;
  2745. //second偏移到第二个元素
  2746. std::advance(second,1);
  2747. //如果两级都能找到才运行继续后续操作,否则,表明没有此路径地图集
  2748. TDOAReaderPathMap::iterator rdm_it = trpm->find(first->second->reader_id);
  2749. if(rdm_it == trpm->end()){
  2750. //表示地图集中不存在第一个分站的路径集,则删除此分站,并继续循环
  2751. pRdm->erase(first);
  2752. first = pRdm->begin();
  2753. continue;
  2754. }
  2755. //统计first分站和后续分站的地图集数量
  2756. int nCount = 0;
  2757. for (;second!= pRdm->end();++second)
  2758. {
  2759. //确认第一个和第二个分站之间是否有地图集
  2760. ReaderPathMap::iterator rpm_it = trpm->find(first->second->reader_id)->second->find(second->second->reader_id);
  2761. if(rpm_it == trpm->find(first->second->reader_id)->second->end()){
  2762. continue;
  2763. }else{
  2764. nCount++;
  2765. break;
  2766. }
  2767. }
  2768. //如果first分站和后续分站都无地图集,则删除first分站,并重置迭代器
  2769. if (nCount == 0)
  2770. {
  2771. pRdm->erase(first);
  2772. first = pRdm->begin();
  2773. }else{
  2774. first = second;
  2775. break;
  2776. }
  2777. }
  2778. //定位数据条数不够2条
  2779. if (pRdm->size() < 2)
  2780. {
  2781. return ALGO_CALC_SOLUTION;
  2782. }
  2783. //先做第一遍筛选,
  2784. //根据前三个分站计算坐标,如果有解,直接返回;如果无解,做第二遍筛选
  2785. bool bNoSolution = true;
  2786. //std::unique_ptr<POS> first_pos = LocateAlgorithm::Pos(pRdm,trpm);
  2787. std::shared_ptr<POS> first_pos = LocateAlgorithm::Pos(pRdm,trpm);
  2788. if (first_pos->posx != INVALID_COORDINATE && first_pos->posy != INVALID_COORDINATE)
  2789. {
  2790. std::shared_ptr<POS> sp = std::make_shared<POS>();
  2791. *sp = *first_pos;
  2792. udm_pos.push_back(sp);
  2793. bNoSolution = false;
  2794. return 0;
  2795. }
  2796. //如果无解且数据条数小于等于2,则返回选择失败
  2797. //主要是筛选两条数据,但无法定位的情况
  2798. if (bNoSolution && pRdm->size() <= 2)
  2799. {
  2800. return ALGO_CALC_NO_SOLUTION_WITH_TWO_DATA;
  2801. }
  2802. //开始执行第二遍筛选
  2803. ReceiveDataMap::iterator first = pRdm->begin();
  2804. ReceiveDataMap::iterator second = first;
  2805. if (bNoSolution)
  2806. {
  2807. //因为first和后两个分站已经通过上步判断无解了
  2808. //如果按3偏移,存在一个bug,
  2809. //即在在如下情况下103,104,301,102,105,106情况下
  2810. //在第一步时,当103和104,103和102的结果被判定不合格而无解
  2811. //此处偏移3,即会导致301被跳过
  2812. //所以改为重新对后续求所有解
  2813. //偏移到第二个元素
  2814. std::advance(second,1);
  2815. }
  2816. //检查偏移后的元素是否存在
  2817. if (second == pRdm->end())
  2818. {
  2819. //偏移后找不到元素,返回错误
  2820. return ALGO_CALC_ONE_DATA;
  2821. }
  2822. //获取第一个元素中的数据
  2823. ReceiveData f1;
  2824. f1.antenna_id = first->second->antenna_id;
  2825. f1.reader_id = first->second->reader_id;
  2826. f1.rec_time_stamp = first->second->rec_time_stamp;
  2827. f1.x = first->second->x;
  2828. f1.y = first->second->y;
  2829. f1.z = first->second->z;
  2830. for(;second != pRdm->end();++second){
  2831. //获取第二个数据中的元素
  2832. ReceiveData f2;
  2833. f2.antenna_id = second->second->antenna_id;
  2834. f2.reader_id = second->second->reader_id;
  2835. f2.rec_time_stamp = second->second->rec_time_stamp;
  2836. f2.x = second->second->x;
  2837. f2.y = second->second->y;
  2838. f2.z = second->second->z;
  2839. bool bExistPath = false;
  2840. //如果第一级都不存在,则继续找下一个元素
  2841. TDOAReaderPathMap::iterator rdm_it = trpm->find(f1.reader_id);
  2842. if(rdm_it == trpm->end()){
  2843. continue;
  2844. }
  2845. //如果两级都能找到才运行继续后续操作,否则,表明没有此路径地图集
  2846. ReaderPathMap::iterator rpm_it = trpm->find(f1.reader_id)->second->find(f2.reader_id);
  2847. if(rpm_it != trpm->find(f1.reader_id)->second->end()){
  2848. bExistPath = true;
  2849. }
  2850. //根据距离的正负,后续判断计算位置取舍时使用
  2851. int nSign = 1;
  2852. long long diffTime = f1.rec_time_stamp - f2.rec_time_stamp;
  2853. //计算位置
  2854. double dist = CFunctions::getDistance(diffTime,CFunctions::TDOA);
  2855. double readers_dist = sqrt(pow(f1.x - f2.x,2) + pow(f1.y - f2.y,2));
  2856. if(fabs(dist) - readers_dist > 0){
  2857. continue;
  2858. }
  2859. //保存解信息
  2860. std::shared_ptr<POS> pos = std::make_shared<POS>();
  2861. if (bExistPath)
  2862. {
  2863. //如果存在地图集
  2864. std::shared_ptr<ReaderPath> pRP = trpm->find(f1.reader_id)->second->find(f2.reader_id)->second;
  2865. //两分站之间的线段个数
  2866. int seg_num = pRP->nRealCalcPoints - 1;
  2867. if(seg_num == 0 || seg_num > 100){
  2868. continue;
  2869. }
  2870. //因为双曲线与分站之间第i条线段或者第j条线段分别有两焦点
  2871. //或者分站之间就一条直线,有两焦点
  2872. double xcross[2] = {INVALID_COORDINATE,INVALID_COORDINATE};
  2873. double ycross[2] = {INVALID_COORDINATE,INVALID_COORDINATE};
  2874. double zcross[2] = {INVALID_COORDINATE,INVALID_COORDINATE};
  2875. int nIdx = 0;
  2876. for(int i = 0;i < seg_num;i++){
  2877. //计算位置坐标,双曲线和线段相交的交点
  2878. std::shared_ptr<SOLUTION> r = LocateAlgorithm::GetPos(pRP,dist,i);
  2879. //无解或解无效
  2880. if(r == nullptr || r->nCount == 0){
  2881. continue;
  2882. }
  2883. int nPosCounts = 0;
  2884. for (int i=0;i<r->nCount;i++)
  2885. {
  2886. xcross[nIdx] = r->x[i];
  2887. ycross[nIdx] = r->y[i];
  2888. nIdx++;
  2889. }
  2890. switch (nIdx)
  2891. {
  2892. case 1:
  2893. pos->posx = xcross[0];
  2894. pos->posy = ycross[0];
  2895. pos->posz = zcross[0];
  2896. udm_pos.push_back(pos);
  2897. break;
  2898. case 2:
  2899. //筛选解:从两解中选出一个
  2900. //解到两焦点之间的距离
  2901. double deltad[2] = {0};
  2902. for(int j = 0; j < 2;j ++){
  2903. double d[2] = {0};
  2904. double dx1 = xcross[j] - pRP->x[0];
  2905. double dy1 = ycross[j] - pRP->y[0];
  2906. d[0] = sqrt(pow(dx1,2) + pow(dy1,2));
  2907. double dx2 = xcross[j] - pRP->x[1];
  2908. double dy2 = ycross[j] - pRP->y[1];
  2909. d[1] = sqrt(pow(dx2,2) + pow(dy2,2));
  2910. deltad[j] = d[0] - d[1];
  2911. }
  2912. int idx = 0;
  2913. //应该和之前计算的dist同方向
  2914. if(dist > 0){
  2915. for(int j = 0;j < 2;j++){
  2916. if(deltad[j] > 0){
  2917. idx = j;
  2918. break;
  2919. }
  2920. }
  2921. }else{
  2922. for(int j = 0;j < 2;j++){
  2923. if(deltad[j] < 0){
  2924. idx = j;
  2925. break;
  2926. }
  2927. }
  2928. }
  2929. pos->nFirstReader = f1.reader_id;
  2930. pos->nSecondReader = f2.reader_id;
  2931. pos->posx = xcross[idx];
  2932. pos->posy = ycross[idx];
  2933. pos->posz = zcross[idx];
  2934. udm_pos.push_back(pos);
  2935. break;
  2936. }
  2937. }
  2938. }
  2939. else
  2940. {
  2941. //如果不存在地图集,需要检测计算结果是否在地图集上,
  2942. //根据距离的正负,后续判断计算位置取舍时使用
  2943. //无地图集的解求坐标方法如下:
  2944. double offset_d = (dist + readers_dist)/2;
  2945. double calc_x = 0.0, calc_y = 0.0;
  2946. //if (f1.x - f2.x < 1E-5)
  2947. if (fabs(f1.x - f2.x) < ZERO_PRECISION)
  2948. {
  2949. //在y轴上
  2950. calc_x = f1.x;
  2951. if (f1.y > f2.y)
  2952. {
  2953. calc_y = f1.y - offset_d;
  2954. }
  2955. else
  2956. {
  2957. calc_y = f1.y + offset_d;
  2958. }
  2959. }//else if (f1.y - f2.y < 1E-5)
  2960. else if (fabs(f1.y - f2.y) < ZERO_PRECISION)
  2961. {
  2962. //在x轴上
  2963. if (f1.x < f2.x)
  2964. {
  2965. calc_x = f1.x + offset_d;
  2966. }
  2967. else
  2968. {
  2969. calc_x = f1.x - offset_d;
  2970. }
  2971. calc_y = f1.y;
  2972. }
  2973. else
  2974. {
  2975. //在有斜率的地方
  2976. double arg = atan((f2.y - f1.y)/(f2.x - f1.x));
  2977. if (f1.x < f2.x && f1.y < f2.y)
  2978. {
  2979. calc_x = f1.x + cos(arg)*offset_d;
  2980. calc_y = f1.y + sin(arg)*offset_d;
  2981. }
  2982. else
  2983. {
  2984. calc_x = f1.x - cos(arg)*offset_d;
  2985. calc_y = f1.y - sin(arg)*offset_d;
  2986. }
  2987. }
  2988. pos->posx = calc_x;
  2989. pos->posy = calc_y;
  2990. pos->nFirstReader = f1.reader_id;
  2991. pos->nSecondReader = f2.reader_id;
  2992. //如果不在,则此解无效,如果存在,则此解有效
  2993. if (IsOnMap(pos,trpm))
  2994. {
  2995. udm_pos.push_back(pos);
  2996. }
  2997. }
  2998. }
  2999. if (udm_pos.size() == 0)
  3000. {
  3001. return ALGO_CALC_SOLUTION;
  3002. }
  3003. return 0;
  3004. }
  3005. int LocateAlgorithm::round(double value)
  3006. {
  3007. return (value > 0.0)?floor(value + 0.5):ceil(value - 0.5);
  3008. }
  3009. bool LocateAlgorithm::CheckPosInValid(POS* pos,ReceiveDataMap* pRdm,double dScale)
  3010. {
  3011. if(dScale <= 0.0){
  3012. return false;
  3013. }
  3014. if(pos->posx == INVALID_COORDINATE || pos->posy == INVALID_COORDINATE){
  3015. return false;
  3016. }
  3017. //判断解的无效性
  3018. int nSize = 0;
  3019. bool bRet = false;
  3020. nSize = pRdm->size();
  3021. double dist[2] = {0.0};
  3022. if(nSize == 2){
  3023. int i = 0;
  3024. for(ReceiveDataMap::iterator it = pRdm->begin();it!=pRdm->end();++it,i++){
  3025. dist[i] = sqrt(pow(pos->posx - it->second->x,2)+pow(pos->posy - it->second->y,2));
  3026. if(fabs(dist[i]) < NEAR_READER){
  3027. //误差在范围内,判断此分站是否属于特殊分站
  3028. if(it->second->special){
  3029. bRet = false;
  3030. }else{
  3031. bRet = true;
  3032. }
  3033. }
  3034. }
  3035. }
  3036. return bRet;
  3037. }
  3038. std::shared_ptr<POS> LocateAlgorithm::MappingToPath( std::shared_ptr<POS> pos, std::shared_ptr<TDOAReaderPathMap> trpm )
  3039. {
  3040. if(trpm->find(pos->nFirstReader) == trpm->end()){
  3041. return nullptr;
  3042. }
  3043. if(trpm->find(pos->nFirstReader)->second->find(pos->nSecondReader)==trpm->find(pos->nFirstReader)->second->end()){
  3044. return nullptr;
  3045. }
  3046. std::shared_ptr<POS> p = std::make_shared<POS>();
  3047. //获得地图集
  3048. std::shared_ptr<ReaderPath> pRP = trpm->find(pos->nFirstReader)->second->find(pos->nSecondReader)->second;
  3049. // 计算两点间斜率
  3050. // 垂直
  3051. if(abs(pRP->x[0] - pRP->x[1]) < 1E-6)
  3052. {
  3053. p->posx = pRP->x[0];
  3054. p->posy = pos->posy;
  3055. }
  3056. // 水平
  3057. else if(abs(pRP->y[0] - pRP->y[1]) < 1E-6)
  3058. {
  3059. p->posx = pos->posx;
  3060. p->posy = pRP->y[0];
  3061. }
  3062. // 有斜率
  3063. else
  3064. {
  3065. // 原理 y = kx + b, y' = (-1/k) x' + b'
  3066. double k1= (pRP->y[1] - pRP->y[0]) / (pRP->x[1] - pRP->x[0]);
  3067. double k2 = -1.0 / k1;
  3068. double b1 = pRP->y[0] - k1 * pRP->x[0];
  3069. double b2 = pos->posy - k2 * pos->posx;
  3070. p->posx = (b1 - b2)/(k2 - k1);
  3071. p->posy = k2 * p->posx + b2;
  3072. }
  3073. return p;
  3074. }
  3075. void HeapSort(_coordinate** pCoordinateArray,int nLen)
  3076. {
  3077. BuildMaxHeap(pCoordinateArray,nLen);
  3078. for(int i = nLen - 1;i > 0;i--){
  3079. SwapElement(pCoordinateArray[0],pCoordinateArray[i]);
  3080. AdjustMaxHeap(pCoordinateArray,0,i-1);
  3081. }
  3082. }
  3083. void BuildMaxHeap(_coordinate** pCoordinateArray,int nLen)
  3084. {
  3085. for(int i = nLen/2 - 1;i > 0;i--){
  3086. AdjustMaxHeap(pCoordinateArray,i,nLen-1);
  3087. }
  3088. }
  3089. void AdjustMaxHeap(_coordinate** pCoordinateArray,int n,int nHeapSize)
  3090. {
  3091. int l = (n+1)*2-1;
  3092. int r = (n+1)*2;
  3093. int max;
  3094. if(l <= nHeapSize && pCoordinateArray[l]->tt > pCoordinateArray[n]->tt){
  3095. max = l;
  3096. }else{
  3097. max = n;
  3098. }
  3099. if(r <= nHeapSize && pCoordinateArray[r]->tt > pCoordinateArray[max]->tt){
  3100. max = r;
  3101. }
  3102. if(max != n){
  3103. SwapElement(pCoordinateArray[n],pCoordinateArray[max]);
  3104. AdjustMaxHeap(pCoordinateArray,max,nHeapSize);
  3105. }
  3106. }
  3107. void SwapElement(_coordinate* a,_coordinate*b)
  3108. {
  3109. _coordinate tmp;
  3110. tmp = *a;
  3111. *a = *b;
  3112. *b = tmp;
  3113. }
  3114. void SelectSort(_coordinate** pCoordinateArray,int nLen)
  3115. {
  3116. for (int i=0; i<nLen; i++)
  3117. {
  3118. int k = i;
  3119. unsigned long long key = pCoordinateArray[i]->tt;
  3120. for (int j=i+1; j<nLen; j++)
  3121. {
  3122. if (pCoordinateArray[j]->tt < key)
  3123. {
  3124. k = j;
  3125. key = pCoordinateArray[j]->tt;
  3126. }
  3127. }
  3128. if (k!=i)
  3129. SwapElement(pCoordinateArray[i], pCoordinateArray[k]);
  3130. }
  3131. }