From deccb290e55154c3bc391c824766a73441df8f1b Mon Sep 17 00:00:00 2001 From: "Dr. Zygmunt L. Szpak" Date: Sun, 17 Jun 2018 23:08:30 +0930 Subject: [PATCH 1/2] WIP Triangulate --- data/teapot.mat | Bin 0 -> 41369 bytes demo/example_fundamental_matrix.jl | 4 +- demo/example_uncalibrated_triangulate.jl | 113 +++++++++++++++++++++ src/MultipleViewGeometry.jl | 20 +++- src/constraints/ModuleConstraints.jl | 7 ++ src/constraints/constraints.jl | 47 +++++++++ src/cost_function/cost_functions.jl | 18 ++-- src/draw/ModuleDraw.jl | 3 +- src/draw/draw.jl | 78 ++++++++++++++ src/estimate/estimate_twoview.jl | 2 +- src/noise/ModuleNoise.jl | 6 ++ src/noise/perturb.jl | 15 +++ src/triangulation/ModuleTriangulation.jl | 9 ++ src/triangulation/triangulate.jl | 56 ++++++++++ src/types/ModuleTypes.jl | 2 +- src/types/types.jl | 10 +- test/perturb_tests.jl | 23 +++++ test/runtests.jl | 3 + test/satisfy_epipolar_constraints_tests.jl | 86 ++++++++++++++++ test/triangulate_tests.jl | 62 +++++++++++ 20 files changed, 547 insertions(+), 17 deletions(-) create mode 100755 data/teapot.mat create mode 100644 demo/example_uncalibrated_triangulate.jl create mode 100644 src/constraints/ModuleConstraints.jl create mode 100644 src/constraints/constraints.jl create mode 100644 src/noise/ModuleNoise.jl create mode 100644 src/noise/perturb.jl create mode 100644 src/triangulation/ModuleTriangulation.jl create mode 100644 src/triangulation/triangulate.jl create mode 100644 test/perturb_tests.jl create mode 100644 test/satisfy_epipolar_constraints_tests.jl create mode 100644 test/triangulate_tests.jl diff --git a/data/teapot.mat b/data/teapot.mat new file mode 100755 index 0000000000000000000000000000000000000000..a0d30de98038a4530ff048dfd289646ae5bd251d GIT binary patch literal 41369 zcma%iX*kqh_;xii;{3rHA`wajlUFM6+QvTkdt?XroEOpK1P@$Ztv zYmxa+-A0onlv?y6#&c^&KK9JCe4p<9@7$!->VNOHP!r#JDp?r>OE7>R#qUE1!v->KS`;#Pv8I3-8r>_SbYmtm7R0>%-%-JKDc$sYtU3$gDeZ${03@I ztN66r=2Kt-?(PiI>@`lPNiHQ?i`IT+n^uDOpEv9~w)yBdZWf*!>GN?M3j8h^7f7mA zwSICo)d5p6CbRUouDp8lzDwQUd!j^f#o?YB-bbdhcFf!9c@^*=MFA;fn@Bt4zx-(* zt9QyCzA!A+C4T)lXYfd!zqGZGTG1DE*vMm)^OJE!Hxby-O6~XS0!Zhh+v)xLCd zvL;*L0F=!+`+h>^ES>PssAEn2)n`@fdVzF1jI3g~$9Y+&TS;A=scoTX%A~S5K2>?Y zc)gofo(Jy+zIaZHoAg+cj;lF}$$cZv6i3`=??%A?I$U0ZK+h`UM~gpj6|Zl6z!jb(99iHxQLw9j0{xPH*<*yo>H;aIiPnA zWl38b^_*@34+blrMz~@->}M;62a?|dsalW1f-c9Q)X|ZUY;vHR!HKbwiYU<$cZ3Ac zFbXfD_AB(?@6dWDviDtewyPE*!niu}A@gU&ok{4ejXW==R+Mm|rPwHT`xNJG1oh3R zFp22A9Q1Ym7~MhfQw^dxLF&LW@JutS)K4RPuK44c!S%X%#Km9I1=kX=CCPYuRPy5sEQ z-l$^l-@javX2p>n@RLY;bF0&c09YDUmdX)Ka?3_gb6C-N4_jJ~N?RAKtE1NNUyg4< zYb|d6(gE5aBlg=DQbuAjhs%t92cFFW$?0pFBv!_EQqSQCdI$Uy8>)zbrCUPT5ayN+ zN0+&=^Quh(K`MUtGKsKxFyDI9ce z+VVTI##P{OnlLRT1rd2$d`d>B{#<}+)w@u{BYZVk^1w_aOKp7_1~Siu_*ii6CCrs8 zFBbFU@v=vwE{T_x>3`SVoju%)X0KW5C$`j0gq9(Ep6Dxj!jri6fXrw0Cn@#Hi8TBx z{K=1yg=Mk$g3ETmUec%bk94A|xz{tbm(k1(*VOYJ4R*&_@wOH!@Dn4M@%xGi^Lw1E zyi*TmU`1+C%(oUArxB@hrw?0aUhklUq>Nf{n4sm8Oh&hpw6*`{L!sXY?cg@WmCBQQ zvZF1h;e%si5%86e!7|e~DbC`wb#*?md6*TuBtm}eMp=r?Y`1Wi|OtH-2aZ^d1N z!^((y#E#a-&OiC}lD)=hW0`FT=7?$AB`wB(dwHVdZv&ID`(60GyHOX`w401x+T)D+ z{4P6i2Z_^AHaCS>-I{^{wjgoVUl-QnTeyVhrm~CGbOBG+rI=3C&b_9q65ydjzlh+% zzA2hJE1&5c-crN#(F~(IPCA_cIdP+#fl3{Y2_D&fX~v_&3?ASM%#27ezq?x`1!-qV zW5R{qd@lmM^`Aq$5N;e0bH(L8u{|85lQ}+m>h_brXQb(ykyQirjApW^G&6y*O?f4^ ztU!-kk~@khOrkG9rQs(*-@lKrJ8;%-`^{+}%Z6K->9|GRQsf8i8K8KhblQP0i;q+G5X#jBqnub--xmcY_y&1^_3{1}4ZC&bFZ%{|_(kRW{X?@H(2=jPgpue-|FB9IWaG!#D2J$w@)E&o?Esct?D`B_*!@lCkn`=s}0{k z4o39r^RRVj)8WgnGG8F*v-ou_iRL`owZ~PziYhjm-3}bL=W$!CZKv1lP}?T8fVT>r zS9F4NE`pB$?~oeU9_>vhxc~1f?{%h4zr_celqNG|34@`X2VXp|Eyurmun1J;p8_C+ zC((+FOc&f^4pGy&#t#&Uu})f%VR93W!1u<3IlyoX;*Vp-&nACF#KSVWn1A^_Z|t5G zHHLgaI@5V}9s~;r=u^0Y!|GmQkw*+u1iNh_cAh0HO%XNQxLwB$p}!HWG~?BV1`A!` z3;Pwc7?B-pJoFo6IId0Y*0}5$eBKK7Kfts_xgZ!KPYjz|#<^oK9(kaw&Pgd<#f;UO z5bQSOBlC|q(jP3i$;GU|DbJLwcaP4DQ$_GrhL0t@d$iV|JHe4lX-_gvgVoM*jD&Zj zt}7QU8pRa%@9%wIZE(|kPqNNNodPuzJ`eCk_rOL7#tKd7;%2!n@$pQfeY|43CgNEe ztaZi34tV%gaS^zDIRq0ubF(F33j&wxS1zqybArp@-0toZB7TWuYQYEG9*0K1gEK+CDy7rmP?m0IF6 z06875#-Ynf;2A)POefaA`^$UZ06a3sWee*+>3fEq>0F}1qv8YSPZ>WQr(k27i zFB?s8DfrgsV#1KP~(4oDv!hv7p&cuXOfEByJ$DadzO&o{c`@615irvX7L1oN=9X-6tIS3@ib5w8nj~E@W5T4f;1;mHFIgp;ZS-Ox!N?9R2KJ-e{aIgKyKz zsS`d^WVRt>JJinlwjmTz2ff4kcu0s~WqZFCB35SxsdovlY}4N&YL+H1v_AwJvx^;h zH#D9Q66rikdAvjNM+u1Wh;;p?N%1jn`!2N_#1gi>Yq)OZXOyb-bYo;E(oww3tDOCM z#G;ih2y$BrK2y%VbN@KTJQk-N{4+vS`miV6v<$^F-@H^kfjHg2&e6D&&vS;=M%$Rq zXkZHmiFn+bhWDbnfK~Oq@8w_)R_+W)WCztvVZtqX609~g;>bC+BCZ2hHPwTP}R zjxoE+F7$%ULEb~+#|!M3N&6Z!!WWh}KYubKAH)`9N;we0^Om_bEpj^0;c$Fh4KNQ# zfM3@h*32FlE|HwPrPi_cu@Vc=!gg*Zj4BTaBdj>)=5H>K7!RPFUHI?zgUo(<*50HN z#KI8Ww~jV|E$)06$LX+?rMc&WNrHLvY#7me_rJrzx#zo%t5C1GUf`o{7tlQnBq_e` zT*{pHI^+rfsJaU&wAthacJi(E55#f#sA<0~c^k+XN_Yt*bbr~wwIK_XqGE!2tDAZ( zS^Moz%D)2k+buj zE0BI#E2DcycMPbE0)K8(ASw>BDj$@0rO|esMVf$|NHRopYv%*ijpsRwmE|KCMO*Y? zZ#5O3aH&&z{r7&dg4R1E?&|V!v|H`V!Q6Ar6dR_q(d`4dw!`j*H>>LDip$mA4bYh5 zZ45<}@++gMM`|pMPa#HszNt4!ZulFGpkVW+Y(;nU%k_JI3Yix3MqO%ks1}d-7bw{M z2zW2!6!YwZEGqS(PlGJU?e6)p`D)uS@2mC_8`7&6k7|fBC)r1j9AUGAEd1X;8NUqq zrGG(H-*i~~t&}Q&gI;RrbM!ppqiE6KCcH^HHX65EDTrl2O zezxus?5k5y_#vQeq9GI4zg_RrbJ!2^LpY-A4{cnwRT*P zn*RLa*S!Bc9nbC3=q0z7mz-VfKdk6ft7RQCo#{4c(VR^k%Rg04C3`pSw1|&s0M*h0 z`$(qQ#IQOj%&WHcUV6Pst)My%K))I6gS7D5qRu*;p9I9-e!XYD;WV%N{6(VHtmb>~ z7*fCfRRhbM@=BF2M7HIu(>L}_i>@TK(9*-raAsd}#hXQjC(0oLpr5>c8vfxA24_m{D7)zN_Qs5^>nG^b=`&PHua>fAB{?r)_XF)9szr=XZ@o(3B^(&j-<^91uoJZu2~r>G&Hm9qcVkX8(1(c;LQz ztQO&8A>IExLqq?0F37cuj|kh*7!jdE&@b*EK&fw;_1N&Ba_)Rq+V;T*8^KGC7Tz|@ z_7k401$*@Yphu60PH?XMTxH3m)u?)=#OZ;n zRujPAAO)c?p9o^ynID*g#mNMGWo>d4{$->uM4MU$o z<=L|6UxNR#=n@gL%PIZM;7`mnq~*fgGi+eAu70aV*6*NIc9Id{@YSV*6bs?g;297v zeW3KufVSTYC7}cR>fhVz`hC(wuB92@HYa&B&n_bPvGmw}%gwH4ksL=ps@uLE5#=fD zTRiGQ%GAo(<_eqXs2w(@;A)lKA~I_%dNrH~VOV8q8|%e(Mb=nDX8og?VT{i*ywbpQ z9`&p9u#DTf)@aH60a0R&U}1QAbu&OgoUUjYPjUNPa=v19l<=^NRd@(|U6`4ej%Kfx zb2>%r{MiqV>+i>gfC@Qu*Qwj05||0mF=wq^P)Sh6hv2vtX$#3TmC$ZPj+f?IqEHq> zx$4d=fRs6K9ofp@q%Cn}lD&t6`S#R7W4r-oW>Tl7hR;PEnMpIwV_jHNR}9&K=E*EZ zO_f~Q3LmAuDb=!`T&oqvo6*qJmY6nvrW*$(T<}6jjodW)&Yt*Q*wC8bVifO%0v@VR zUTlMa{rL)g%=N@Q{spP6ia3q~$xgM3-KZM#egEeBE}xgs7M1?Rq*ASEvAu4PTzh7z zo6t%2v$;M#N|w%5ywKcXp8cG-2vryMRC~?YjgP~n@y?Iyi_z;-uBu6BeIE+)O6F;e zMRvjPx2{72z4-T_OI%u&FPHC*L2sYed{c&u$joK@Y*Q?w7QH{us8?SA0_(31Of!FJ zf6V!5dov3t7lyZToQ=#VphLe-;84_6@ij}?jdQHZi3<{d|JsONgOFdS=?Z-iZL{In zSmAFI-J@Ca>gpiW!*qfluc7m!Grm;{@?YK}z!(-JQP(5fE|0YYr!D22#B6#eE#~x1 zH*q5gPp`5+@b#4Apz6tniyVHNU{1!R(Gbx0^5fSq#S3Gnf>a-bvTASz&-E~^eMAIv zl71}jo(=YQp6rTGve;pHV^gDZG5;jnsHw5#*%nHN7xx8O9geVs54bc+Fw!>3Vu+%W zrTUjS@h}VIJ6<5Z?e=fr>!hy&W#^8|Wpc_QhsMwq#=5EU4-R>k&SQafaK8%&+5Y`U z6+Z4OLdyAAcG8ZL*tA{l>#PldjPs4l`mOyVu}5d48KRKp}!B z0oLOME)43l7cO|R&)$r@eClYR0J2SFMbuqzqm^?Ge-1c-cQ*^}XR* zEOuvd0_SkY{j=wo6f3v?+EL7kHYv;e&#Ei81Y12+%|b(Q1lKWZVuDj<3$i{U>_pcd z=J;B7QNM6cXB^i^+AlI0<{^d3!v*G`JdCVFtMfrHn}GA|;eBx4ZgAtYW&`35bRu*F zG4PUni9*=4mn^1^Z&ZkDEdZKXIwMgg!8s(Yqws@*!yixT)KwrO_FyhclO8{kxxyL?3 z?XgZa;q-N(%WCSZqS`0-z!tNXZ5{=zl~Ulxj@<%>Z62DMw`rM=b)|f9G?XH@O7jmh zgeE@I5L)K~Dri5)t}O+YnDH*A8H*l;I9lYyFa^PC-&N5dNDSlE!oA;Yd2=Di$2bRS zpVTHy?gF#D`Md-Wk7b1vBc9Fi(?#R?NT|<1d>^aKsImFZh)5wv0+s=fZt}BElQ6nA z{qQeqqw&~?WJBc7ocrC z#MXAT(so*2wIHJ0V4;9xX%pesP+Xhe&Xnm$PeG$CGzf{mt>akUlKCn2WIq=>jFLyL zvgsq}wLIA@ZgWAjD z8$IzamE~wI%-LrOG|#_}LX%(92mRFU2E8(Pd8?gGv-!1)^^YpcOIJAlL*~%*0A8x1 z2mcSFB}JR5-!Ntut(^#7$MUxM$r3D{VfZ8~%2;Dw%YLE_FwH^i;)z}BEW0kfE_25V zF+f0?SvI?ZBoqYsNHCa09N9HcJkPVMJp{<3r}m1|p(G=XkE}=*Qt~Heb%uvrw2!*L z@!e5<6q`ZNgG!Zz=O_Op*as|auE8v*{QKJ-CP7I!-y3U^hH;vrbiV+|$@Y(rECfyl z^`1JPgEKNx;X$ePXQy}3B~ zyKJ2*Mf+C41C1Bh<6TfscItDlFtz~dI#DR6*WuJ4ju%GBxxQK*E? zJ+99VH-}Adc^5zmPHbeejlbhwjEL_spDMl-r|Bq<$6_!M0IA^496g2DfPG2G|> zkKy(TlP)Uhb0A9!kKa&nekf?5ld6|1#iJ}EXl(piDtqCAsf1M6WhLRt%3l#C&y(KS z+vk&aPmB8e?D#c!O=XMB;n&$LAJV%|P$ ztSdeI^1QY2!+LY}+WF(d?@3+B=SfxiIsd)>q-FhxA>9&SUO~DWZ`(S&weGw-ub-iR=rAhRPdI0?z241QQu+piD9j1=8AzJ`Wj(6bD|%5aJ(*om z6_+r*WYL% z&K-3BqM&~7tqK30b9v0W+Bb9BWyI{wdtPQEk^3B;wR-yso2$@(Rb!E{AHO2$uJ4Bj2cnR=q41~!JS4QrxPw~eo2QKAFUtg?2qXHD+gt7ufJ=fgx*#w z3O!Hzz0ii~z0Aj324-F!+ds|G0s#cFED4N=N9J3~kV3Tu#D1~UH_-m2t+bPcZLD)> zXfpsW6T84RUG0})_(5X>9qOwVGN0pH%pyd!ILWWyQ9U(RnID`Y*q@7~ie3CxB#BoB z!-;6Wd3%qY@RO3w28}nIAu0U7I5BNG89|oSdDEi7v&aES7m`G|QBH%t<=PX@?eKC3 z4juD3=9^4aSZZ2Wi9eG9i9xKJ%5Znq1@PSoL5?XachNX~E>xdAEYS5}eea3r-?LBF z#x*)b-!ng6DHV7JazVW%b^N7ks#9MN*}mT+T&qkna30?g+WS`zbJ!}ZjkUb7h+8pT zX1$n|^eQjSg)A?9zoY^e-VP5O9kFW@w*Y&U>FO2(-jwe&CJ0I7#@T=@HhW**V6=Yo zquKukHl-=5gcqadJky1>nlekzDTK2H0`*VX9^BStKwa+$kb4!=ML&2Ip66o^rd%Y4 z*pN>}yzp+M7DXVuzC8Lxax;`ga0UI`$=PI1B_4;eWIbSV9&8Rh*8b5p_(Z6X9^f1` z`a{{zuZ>&&dlJR5#{AWs!16K3Q1uSCG~<$^Hh75-Ydn@eI)TzzJC99hr*{~AgvJ+w z+Yr6-{z@;#r+ETKIufPHj^0$Ip;Ys9+`XUyA=hoa1H!AHI*OxIeujLer~%oCTo~mg zQ_$$nZ36$#ZRUc8b(2;E3-cp8fw+MMZew;{rY4HwJBW=Efm@rd>b6Gic3w% zQu8vuxI7-@+)Z%FA5EM#STIYY%8K{Q>cMCv*iX zF@K4a*Cjf}(kEy|Kys&Dm$`hJ45HCf4ez_zliU$Nz*|rVe9=6N`(R3x>j%J5&^*>< z*O5%g`}!{PH0=y19sMQIP7kAwSr}_jm%BjGG4{0MWRf4Gqd`=Wny4oh_lTjawWoXy z{LduLckez)a6!Iz_%*s2?^P(!XOu_}sgIBKT!SV8TTntg2NZU&&H>J>@NxNv$u7wT&ilM$>SJR-OB_c`$A2Rp2ULE4qC+_vaPckpwc5Iu`J#wdbe*A3Q~W-W7-(GnEJ zeM;mOwluJXeNeIzTw*+3mQ3v2n+gouxyYs78<|(YR)BL-2U}$3y}H~|AsLAa?2Z9a zQ$e@AkHjFh7|)>#SmZRTLCDIrAgvrs6Li}g+02$=|0W(?Wgi1gh+^jbmt{)E???cH z&i_%AAfnOqFWaW;6kM2$jH~M>*BppbcGmNRYC)szmv$p(gN!JS$?D2~spcF% zL|?ah-OJeCxldGrh%B6B-4qU88oh^3W*x^0h;NY1cO(({Jf_#c>%*7!T)I`(kbmU% zn23Kkb4CxWa6OFKMMBN3IoP63?h?7xS+bC*Xo__Lw^&p%-!8X&39!Z`cjkhJNjF44 zp{x=7P~5Bw+#h2T> zsgN-%my^t|XpnQC&~)K7Os(<}a!;CDI4ad}Vb?oYbSoue<~TE7;Xx&^8zeG zlnv7gx~FZm+xSx%Pt*%RoGAYpMsIdG628LW3_Ll^f&It`Fg{pdUHiezAaV@UWB%~3 z{SGc$d^~>ztajYPq=Wm7x7@OlFHWC*FA>@x=67V-#0Xr$w)Yl3PXDI~7_)mQrk25H zJuEB;GPK)?Q2{(m@_S7o{1lBZ>fwEM7qJkD`EKTQ{@};Os-Kie2gg`dSP|-r&!iOi znnlS+ge%k7VOQP+S1yDcdbjBTcMs-8B?M~N9?sBJ#{M$ZrgKEpP%)^!gl6ynw(;S- zT266#o7$gECrV0r@bSI$8OVQbjnaF2_|h8G2IcbIwb0=Kq=_sou0au>l9Lc=GqsG! z2X8XQ51VlolH8%Dhth+rk?o`N(?Q~f2h~#aQ8mUnUu;-Z3&x`4;zRBpK4$-iNB2a^ z$@6yu?;xnVw_R-5TYSrhc?W79w>Hc$Zd&ez={VI+>`6D@g!ab$@BH`S&%!j7o5qt_ zdOpGT7r7t{|8eX+p9USE7WMb)Z36!O;?_(x@LFU)bdn<`PD&z9tPWc_1$E0e`^(8) z_J48xfU`QULi(86hHB%-sQBWJq=-x&2}*aW_?BrQp|9lcv2i^N2m2F!x}(6(vf&ao zI!*<{FjJe&bW{1o zHaPyMeKa!-))@lnt3J9m<9bIO_mul|U@ahuFTz_r2pYoFroAmOty2gSap7sMgp}O@ zF=a}E1@_`R(}x+J@)tcO54fUGf0++hW0!Ylvk_+&i>*M{bCFp_Iitp*$XrTMnt9H% zH5@0si{%g>cNG4#LFT@1fu2Yz$u@v@zB0;$C&rODIMs_N2tnbZl zmNe(&M1liycCLg$NcX5EaRsQeCTKlt>l2$7p3UgrOZ6qFpBou|$TsuS*Gq|XQ`-(f- z{lKVCJ;`lzzlT)6uyq7bN{NC{;rak%w|VxTFmL-)k;1L;i$zFJtox$H;3VfZq&{;4 zK`rnK`RI59)e8nI%*mO@^XHiR^%911ML(JA4?z>C!EDt%6^(|LxO4>K+Fqua@?~#z zjFY%XsEmML1fsw}m!*ij=G9I!to{L=gnT*WIGFMVG?y9bzfRBG z9<)5T#6&aAPs5*@`|2_~H5C^$G^iLN?4OU4&5p3}58@ub3@VQpv3azun{-A#LH@H~ zOvH=dVn0*=9 z@1Zoaw|4Mw71@jL|66;Q?uFgc!slH((^ISxA*+}vzNz7)XtP0{{60o+_d;XiARml6 z1>*V}CBS7E7Y+r7iFGSQH!d3H4SdW}mS$U?QHs7P^a5uQCZdTleSMlbl5&M+g&v&=H6a#OmBg~t`CUU8n56bpy z$N6~zO5Sf38AL1w`E)Ea`_|Ex3gwC|si#6p-|B(ocMj-dm=Bij!2@~zZ+KWP>lJqu z!Q8c#FDu^gGR^+IGMH~qbhGdn;T9;F=&U7h5N)O3Z|jQBorh!sO@`yvbcN>H zzg)ByaHJ~#)cTvviQOOOr&F1~nSFvF84H(va|7Y6t)RBFF8n&^Y<<;35xZcu1~)$X z9@MPWtuFAX8Uyu7pD|ox<#D@T{@M;FO;=F|f3h{e!o^oj>!jy9GWS(RK%@Uhg#X6d zwUi}>coLv=IkC_h=`em`{WelJSlu_BP=VK}I4MnnY|+t-JW2F!jV z->hT3Q$1MRSm7?IzJkU|?wxE13mna)NFr}2Fo(9i=W_PG4XKKf-tl~anY*R7Lsx-! z>@9(8Xu{p)M~zv>{r};t6@4R)<(t*#!&!Nc^4YT02F>uEjlh?<6-Jvx;homfX+fY8 zTQl4DYxZYDDHe$(lLQ+jMs_J`#w9U&;N)R(V3T8hxUtYkX0N3FSrtfhX*4dPlr9=E z`=vWJ(W{DZ^)Sytc>?ukXYTkx<}4!`9%-OeqDuI7qVdHxuN~+sAmAU>DTfO6+F9q@ z_Rj5NoG5MuH6cfiw>mMsCuR4Bhg(@=)5X6-JFx8uVQloL^~-=o_7LO6$B@A1rW&-n z@B-F1XT`q1Fs?;3|vW-ul7xQ+wJ~maP)L4eC_Htt0H$S7qJ8p%x~kjtkJ1FG_2<1~9vEi3(FR-pQZ%~NRr;y~m+)*3)@=uANCz51~ z#wWx;*9#P|p9%?;hB$`U?X2ek6!&!HvJQO}TkGGySk_9#6`eQ?nSH6i=kMI+fzV*W zFzwo1I8TIaRsI>ob7op%JNQq6#RyzTqHK;<+ho8CAICTyWse@>z8|A*R_)wh8=ZNY z5`Q(Ri93L_Dw!W-+*=B*WfU&|eBt6-F692X>X%`5M!zre#O#bIhkyHR5s*{fJ`r4& zr2>Sn))f$R*JF$ z5*G9Jx9CE$8wPw-8mjH3%XU&c^D3hYeR6T`8KJ){d&>a+4XhUHvkzO*PGXb?}%+Zx*SBPBD}(7KR!U-1m6OE zj<9rDeSQKMHg_ON(2N9SF!igE_6SbI4a|ZB!aMN4Bc5H%*(6$=TwNsSE{b}t0R1pw z{GUpLwq5|eQApr+EoyrHko0a~g;sP>#0{tg8tMh@c42rihf$5+*}=8|>UYDIHwJ${ zXk8EqSJb&mE}?7wc?@qT;+CeSvF{uWzQj$sXMJm9M%eJ+kQmbVhctS@?}Jg)#TJ5S z95wuRZt|Kbh7vml7!}4eJmoz{cgZ~dAo!Lc>=?$3FE;NL`xDm)*k$afpP&IKhDM=S zJEPE|ORR)T^Z|pSB+&y#+s@jfkq4Y>_I(@*XkKQ3s{w1F@5XUd*KXFGYF(DfFU>t3 z?^wi_c-yQApcy`5PzxE~%=l-7Q@YGV*B(mmf_(HVfA|%Kt0OaKJ`SNlmv-^EOQc~< z#;q7oqbFal>cvq`&u-f0^FJ56d?X812!pg&XN7qjN+vrw77)MLHwb~*qJ=_^S%%K{ zlr^`#`|L+gbd64ncyvf+yoIkTiBcXjvt#oWdBF>+P@%I_yvixV!LL;|DYm|XxUpK|*U^b!sI z$kDwpUL~HbAlwiOhWb6aO0sTb3RCLN0zzu^WIA3=p&3Vj3W%)TiQmqzq@ z(aW@V#l+*R1wPgZ!|N%K#8F<7xRlO(kQ2B?Mk7m9@Tj%6dB1;xB;LTGX{XyjL$g>- zc_53}idpM2w z+Hhl%Fin{`adHuu>UCN3^U1Y6LS}EBiGh4oGakviE)%-+n(;&Q3@_{5|6}kzt+I=t`=oV~Q#3!!`sJ-R9S3J+@P`)H zH^x`EC=h?v>weN^JZWTq0GdfCue{gw=*{riUdU_q$oWQ2_E2f%5kM2O=W+OIb*duDcgdT3UVZYCJsR~?%YA+P_7RYuJk?HbC8DhtT7)mf%OQntmxH_ zRAFVRg#=GzN8XM?Ppg0KGM~8UpQNF9nFdn$wZCiT7OkP2le-3e&CWM5?)4e^)q-4R zCBeG%F)A%-K4mN;*L}f0m^idWG+^j@x@;#0+wI#AfJfFnF#>VPNOxX{{dNR8ypA1w zjlg$Bj`?ZrojSM$yZsifLGHdyyox#rY@r?&h94*fxQXS$)FW7tq}prl=_=46X|@7! z-GFME1w;@k3zhX-;Plc3VQNs-tE~_Kfw=h z{5{N}t{T|0JKgdDPsBefoeo|X8?7EChiiYGgcPl-(Z(1Yuxbtklz{7fxV@miBX$cZ z0ejwS=A|rz&B*La7e}mpsRSP4&}z~aSG#t@y|H;JbITqN8Gg!Od&u+!>^s>$!+ge+ zlz35D7Y|0|%C)u;{rV0s4sHhZ%z2Ab%sf1GAgzj826IDJiPAjOe+d*)N3ir>_~U9e z=5^x=>wdI(_JtEW6yyZ~Xqa-p)+^A>Yq0vTs)-%?rLDKC>c?v;XVAZv~oNr95T8ud;;}$Tv3emV#J(J3>ZPYb;a=S1Jucn z;ni|KD67d2B_OP8w7XXj1)@tDcT4|l+NO9iJHN;e7|iB#LS8=A+lOt4F1zeIdyT+U z<)RmXH~!im={ua1jB-FV16gvtZ6><9XphDxeDv+tux1i@t$o_4EiaqMd44k?px$Y^9@~t0K7$aHBDH7J zCaHu=Mu5`b>W_%qyZ|I)#*ofllefbV@b5u{qkbl*5pV7dB1CIG&~%8Mi2vt>jIi*t z;y1VqRk^`?L2oxq|4Q2_sq7Yuva;;#{8!pM^w=|R<-}g;x*sjL{4QTgVEpY~%+-md z_OIv2FZu;lzRLY_ijS4PA(-R^9QOm%-2LdCFt7#;*PEt>4#LqNe0seM6OwoY%e=~S zoTL^$xaoacZRx%A+7hQkAMjJoQ0P)byFt7h`@yGM24J)0wA zjpS;n_?36%SFO9PYd;Ynfc)6^)%$C`b3%-VZ+BjK*$62`cF6#Ajb-K7#3!o2FxWHI z`GMM;OJ;a*g6ltK^Ka}uJxp9hcAuZBSM?d%S8aRw^^L+1|9Jvg7GuwoQ4o}JQ{r)5 zDo$cBeOxktw=T88=fs$iyXD)&hoV|2aAO*T;#VpXAIz(1_G0Ky!zHIng!ix0lK*)E z8XjgLO3eCR7ZWah(mUA4jvBuQ$YaO6g4JQdwIJ2*OZYE`s~lkvH@c?;E{`K+VTP$? z8L}g)bC7q=D(yTx?E{|&Yf8~RCMas=jt*K0b+*l~JtgJk&{cfA%`YnxKY|z-^DpIsTI8o$SyYHopD(SMSmLVL>LwIxym49nkhk;jpL4f0SOTG=T$q933R+RSD37lg`o2o#WH& zAB^W+Y=QkeCCWKqq<794NY+`^MHu)#Zt`>V3Bb6NREQ(mnPrM3Z%@qTKRIh_L?A$I z!!Wk`r%)!UFVsd2(Epb@axRb!) zirILu2IhN{<~_P=rsG2!9LXw|o{dJ_8WR3u4IQG0+~OqXji(}S+qF>)rSSE1@_+l* zll#H4BAb8e?jG!~Yqgnc4F>R1PS6W!6ce`aub(;yKYRa~{-AGm55@_a_+9Z0{$)kg z?1=^FlNo`B*J`XNCvTMV(zs}f7fYbPB$2eq)Ax&91>V@b# z$#l~d8INyjqI)FnAY9w0etrJMz><*%9kf6om&IPPnUF6s4`WDQBBkmrEAGcP1$o@u zc~EuT+3CS}Lp=PpU(Ztt>R-&#;2R%hu{6GeHOP`;^+2`0+xBnTa@B}H8t4DbPU{^B zQD;BPCm)2>SEN|o!^il3;NAfi^JuRI?j)$-k=mdCT8iY{jdBPI4-%{IBb=dq*B*hh zk#N$Ma$Fg#gt5EYu&4c6ViLmY`iJI`^Q_&qi{+rJL8p-QV2?5o58|z&&vI0ECuRRe zGshR?o~dT1!Dfn|Bwejs8?LHJI0}lUJahiU9VAvVH#nI-4uI0{k_q-F>9!|R{Db6b z9y9%dt=mqnYhTAIh9Ta9-n74(Rd*;0q!x13LdLqT(Cp~d*FT!JRrjA0Lwc;6{YT** zWc<=-Z&DzY4{CC+Zi{)>Hou-o^Q>AIBOczDJS_5_jH;0!_(lm!0CVzOw4rQRGd97^ zah)2K=fT*n#R+*nzW$SOUU}Mn_NTw*&U-ltvf??)6_%vD@ZY78p|zu5y>t4u?|Wr3 z=i9e$Dfh>w!&ZF~6_1$&_ZpM^4Oar~`aDPn7Cej~W;A?7@eNNKpW?l;jJrIg|z;3aSLrM zV~apz3E3?_%u80TP1G90!=}c+B0V78u4tvzhWEfbXE7{F*-zWB9AhH5QKsSJE#)4+ zMf@}$SV9qBOiKNFj2W|gPcS*@SR-z_%fkvNKQqa8ccCB`zn!5E2FsYx)9WC@`|2Pu z*v(BjPsq>>$AQWIKhdz~dvSh40&uX#2z0C7&Xe5@6CXe?Yl4s9g-26nP<3OV8dGkO zj5mqTw{R7w788HZy?s4op*?i-H-Im?`hxu=IK=0n{zuZKs663nm7Z{rSvmLNEn!$O zh1373t)#=FI?Z3#qEMFin9oiE1m@SyQNA0DkjX2`oQSPyhK=G#yH+fzHnW0O4V(I9Qh?OW-7lF2qPKuh4v2Kll>uW7}Jn%oO^&c*=IRa7U zDt0g@XIAnRe0Uj)I44m{sumpwq@ry+I6v@fj{kWVI~w* zrm^8>%5q-=JXYZd>$9cfQ+AvnSDgurYca(wIDV=#+kK9I8m@VU!D1F5sw;$INNh{S z=#EE_59Q!Jk&9F-s`>?=cNuLAjb&za8_?HJ(iEQ4ES1&VgUUKCQSO6$GcvKIJVXtR zpiXYj1B5;!EdddO3AHin-4A&R+5l^Wy8Prlr9RqlSP+EXj0#u0B}9#^o#ib1S)^}y zaPsGL!to!Sw*LQzrY~`1`howh6iFpXD#sF%TnV}DE22m%Dmhn(5DU3yOYYo+B1S@p zIdZeP?>qN3w>dY*?B3_`dw$RJFYL3=`}2Oi-ml{?vPvs<6I&xg0F$mfS!xL-`vCHc zi$poDWc~zFy6|<`cmDk3l|4>s_9ICsV6TG1_J+2jldsN?$DZG|8&0YhxdiR)zruyF zM99j6d;F)eStHJ*P~lcu;@gwgGp^VM@t|-ckduHDQ`imrj&yY~?-g(rtED~L7*o-% z6&fIDg=Ma_^!nD!*&M0Odrvu`)CK0}qS(^%+?Y?{w?gLHnOC@zpIZQ3Y3;AA31Y{= zU&w`us{OkJ?>?0n-VNkxQrjkE<<>yUBGd8}TfjBXaW^=UbfOJgUYrfRp`UT>7jQp-E=%~iC%u+D!&9tW7owRZs^EK-{B->?vd%u!H#!(siC@`jdYv+99c;qo{(b~T{QLtCdHyT`&y z6mc8n2EJo2^$T3+3hUVY*KYJ1^?>o)c(;slSl zP(h=mXm!TZ>IV$Nj79%!MB2ieC$L}Iz8h_{dsZHgYK$W+DstqU=Q*W$3XpJi|fRw@VAnEml{VqhD)Vm9W zfz!MC*v;ecr{8fCNf$qJN73Dg>i2@Po?P6m*nRRwZ5fS+7Cxvujfe}dZ_)Z?ww47N z!bZpSd>hXY)?vIoMR#Nr`TU*3{ulc*z)(nI0XKDH`pGLO1Bw4B^MTFTYHHzu zl{tMsT`C$nKAF((Uc8{;cK!F0&zIs~`SNzZsw5g4RPe|FI$6De7JZ|nXapkmJy-{YiZ zqN2q|r?1s|VAD2R2>KM{pb9J+am`}%?rDajqZsn4_@?Of6c{ANqm?xubXaY=R0N*~ z@qo&oq}L!*%t%EqS!|Nk5T$oF(f>@ikL;k~g-LENqtH`)>YYU{jxjWA>(MXE9HMU& zx3sOiZXd%1P$=L{fb0KV)qKv6UjvzoJGut_Ot~z<=kiDU>1VqelAH+-e+OUCS&V!o zns(i;EaH7|brQbqm0^X5q-$Z@NUCvbagNe(az}AdQKjX+I(3Gi*|k4IT}pH4r#Gxw z_1A$OCHy@Uvh{qF5H%5V#3{ih?@DyQu|7L~7|Ju4N?&6w?4{aVj7Or&f_dUMML79! znhhrU>JCrk2hnw?P9Lno!Z$Z&rO2)F12uR00(^Y?+6hL*`AC-q3q}Qz!<-{n%Q3V&olxYf+ww?`EU*b*EziyqR zI60JIE0ys&TPT5>+T6wvKP`dgK=iBPlWZ8z4#Cf;2AS#La;v0_Rg1kL8<&M{bg{^4 zEym{Y#pNpS_}FFFYGJug2y=|^Y1U8j1tvar)2tEK%d{$q^|lnA1Rm!-i<8^Lq}3GP z^t^SQzP9XM2OA6g2jR+q93pq-&cP)>rxnMO_E48EGRLcy%~v>(rk8o()ZM73TLd+> zL$p<)*kPfx3J2XG`sI0c%M!@?#8+uz#r89PB3K5dwg5&iaD3|UUhoo|z`plA8E{4C zF1*R%F#1$5OHbm7@}BN(ie)+Y>mUKtiN*CDDC1Je=L&E)3niB7LUisPz;%c z4$dl`W|SUN4Z)+o|N0M|xP8HNxj@bJq;)Aq&6TJ>0_wSTbc}94`FvU?T5hYxq~LuJ zI(bV-?_fzt_1YekGvbqQM+h%8h^8YI)EM~l$ny0&FBM46MS1)&tb41;5uP0jiPSU(OAkCrpUhM zQ2#gmk&iKco^{cAZ-X)U%JJ>}Ov`HS2Ys1nlQ}7{`33-Drt#xK2@4G@fF05!tN#RQ};b>bf(`2QzUKzgnS%VIxIFd+UYNu($_&b04cu!oIA z_Wjeqck04m+HA3Iz1;NCwbv8OiL{HcjmQ2iDI7u8dP|m!61aeyG*g8<5k$9b)dE}$ zCx50Xrro{AZ@R&k#%)9UTD0!7jJL@9gn_bWfFlPwf%Z(x z(^cz$u*%LXFS}EK{l*q&Dm-;z7qSeYA^58}u>IN;oW!$V2>F?p`?XV8)vbr_2YH9n zuTQ?qv5iv{M8F;w3$Uq5x4dH>YJd)Jm`X2EpM`_B#)Kytj5Simk)nvxVhnLD zJ@{++f3LwxdutuwUBAqgioi^daxG0<;YwV`gjM|0G=@P~ZJZ&@pL^Gc{%c*V8A!fN zvU(4V< zg+aBSkUb=KZ_E-wUVCZZKXEY|c)$Xll7$9Q1?OY|OAeMjgCekTx-CDURbsMCI&8m< zt|Rj^Y78XBZLxF?$-1JyVli{|5cd3X$mfAc|F8=b0o-Ce=_|WT7K-B%s>r-^-+%rX zccZedF%uPY3@SYtLI0i!S^{myE0ZApO|D6c*PX!Ax&NWZ@jpu+Fzrzy(5r40h_ywlA$I6zY8>P$rYO2yo<)SCh7$ms<`NqhgB)zc&*X~ zcYxN`mx`|<*ef8V4tKm22p>TW9J^hvwLqvV*d}{#EmN7=+HI4uI3HD zT81#UGIune=m;h#d1fv&y@d{`3ac zQ`W9OqPI!~bQd%5ACUc|Yr%JM(kB)BCjdzj%+vz_eaEAj{QlVT*ZzvE-TYnSt)660 zJ)weo3!Ynsln`h<=o^AS2XXyk6R4e!ajlTCOrcdT_lfaNSnqJ;D!RyP z5bMS8r>II5K(-7G0@&+Dy`tLTIzJ%zCiFaGrBXDDqv8$UEOL*hvk+%?LQStscP(Yj z^jzRfQG@y?Ou(Idtev-9r^}4hn<3`~=-ygDdOW}SIa*isH|>MzO~EY#a`=e=$F?`d zpsrGI^W5@p_YtzvhAsKhvgQjy*S=-YdmAvV`dplTpT7{Pa02VBEsux(TD9dkhC$n& z9U1`EfdlrTd_lR3)Q;1%+fj|8hCa9*+5*Wvxr&%0E(G(UxKxQCFZj);vD|?IoW% zmCCNL&T$)s)idoOcMqUSw8&6IiK#EkS$lBl5EKqLf>nr^v6>eoIO}>}Y)`D5;yARG zqjx2044KG!noqvBfruDg+1F`ed5)jaLS0z}K?RV*`=^&pM2Q+KiGDyaOzpSnNU&6i zVJToEk(-Tu4uY}Q>oI3!*nM?2HLnKbj<^b?V`U~6f;{U4^%;B9_Aw>y!Aa_3_8CUopkW8zSxRy9GKQUmuP z|Ej6Mg2N9iQP&lf_V-pZ;4b7DuL_nQPl`2(2-qmLE&*38^PsU1&l|fYb$%e0Z4uOL zz3MNvcgiL`ubGgP7WR%wOTpckHytOpXAnB!25Np0UY^gdkmX7- zu({-3?5}EfKkzc;>l%4*w;pTxY%paDVBJT6<9Io+gC65N%(A6^$sapLBO>pkmkke6 z%7qYl5}M1Now<$jCxw56c)XQ*JEQPr&GzHzC9_peheoEdFB zbkN}L)rVU;vy~@5C}56*#nQ*XiNiIPP~mpV7qal!Ik^Y*-1M=6FPYw@NX68@y7j}U zn7ONd@Ttgta?W-=EWIRd4e7@6U%W*SC|EPxvb(vC%i$pT1SnK;pml<&pN%U)j4~2S z_ESoF_I~-lXN)<@Q^RHcQeehf3rsCQ{l{a#(oSug+@o zH95bu>6ZcX&-(vllEbo@uTciL&=Ol2a<6U~fJO#EZOyKt6dF>6AQX#f3~0F@#;zKz z8`8?%rQOs}Li0UVg>*w`4W42ZrPqG76bQ&~J2^P+N8e}Q+$3wDNw6Fuy)XthFKoUg zH&KyYk0SnQp#zIVpbLQhU!CJosE8%ROaLo4S*E+2T5A2^+m>NM&a`eqMcI?(-(pck zj}ylV#Vj898u1iI{LyKzpJeAZu_f*{UEi$lxM;tskSo6w-E`oE%>!ZA?h>d=-e9r@ z)}wzPK2Q#Ql%P1>3^{^iUp?sAmzKWs{+;o88FM$b8>Cs+?h{n})(aZ7w(|Rqvzoed zicNj_<(LH<%=O_^*5L=J(IJkQjx9yk$~ga&AHZR&yHksQpx7{WRhl9`2u=GCB`ar| zTq(53rJ)4#RZesb?0Le9n$z{)lAxLKZ#dGciUHE;D?r&9rHuZX^nvYA<(Fqn9}al; znmcXv)Nkvvfl^|4?3gY1x3wxrv3~4JTgfJvbH&H4X}Qu#?k<8W8>z&WH~x482Cg-^cI+Ly3Z>Z&}l4Mf;-2p)rd%um(!c91VJgpwZxr!mPYOvAqzfy{FZ5hVU<5)s7w5Ug#=za}oQ56Ae2E2TqKV*Q8M{xKUn7AhU>m}*5F zs|s9A*o2iIPw-qw>tM3SN{Oh~8}GC6Xm1gP)uivTyN!~)(VY^bW)b@npDSgwl{X!t z>t@|{G{&DbqHUBHAmE9&xBCNdzG8R>WgZ6U^3E#>p}N6Kd-Uy_jQt$)s*P3lPsj&M z9~=LDP>EvK1az8oU2woxA%PQU$kjUzrB9cBs}j)OF)i?XtwvPGs%IH!L7i*O{qJq* z?@W5!%BZ)$WCm55a6HA*hvvxM*bGRGX~d9C4^BSu#z;^)rMPJV1gBmxeE#wkptt2q zzt}zgwaRJ*oX=jgLMIcyuC(;gS&Mybi0Vq)1`7&NZj znTM4#rz}ktKiZ-44#~dm!6d4A+qn&-+WSZHc!R`Q-Y8#di347D#0Wz&4rZ}^J~NHE zSubzf2!t|;QhOO}BM5!*7R*<4hGz@h>A}UUfvxh32e*M-A;$H&R9Wo9(Aff|9*3a0 z+c}pPU|1LZrF5kl^~-?s_QNF|dW~zlhQJV?+x`4|(KB)u@yCGz)pMj1Pp(F%T!4?A z?2}(5VN%@WHZl^A+aZPH&%ue!@Np+~#J(m3P14-}{w2 z0(`9HhL30t%%gk)$RWXiQh&AbwyxU^O3$R>L1dfaNu0sQS2va!yRM3pNQdx2XqqKH z|L|Q>ic-42e(l1H_8g(rRDB3`YSp2=qKUS)&K1t6=lYO>c()oI(7t?uN-{e?$UnCI z*m>TvHlkJ8s^;3s+NDLUM)K=>Ltv)&DeC7{68-fcQqyQ#ZbgdNL#VuV0v^wL9{k7K9*5F-lv{-IEe*Ey~r5C(Z$b@+r_6LHGl@=xw42m!J zVmNJs6$cy#K{`O=niA+IxO=<5t1+iXK1(0VFkmkE`7g-bO!3oU3Ho0H*-Ud$_3duN zH1L6LkpGkOu~@rKd{c|8IXuaO^A|eL4)Q(48vA%?ED+9w3JN)BROAi(MxY@&)*buS z%~)qX(bLF{WZ+DIuko5K@tIa!J7MtH;ES~14?7+9Q0|tW!zEUzCV@yw;Hzeok<_xL zabD}IKUHRW2PW8x+auphH;Qy-lB6{1cp29!j2e&3or`nIKSU_yovtPNTI4B>dES@N zdss;!pOR;DC-N9io583$PFn9+rMBR3hh^TL%+%KcLs)37Ry4&S$1RZPuosrzh%5j# z;OgNE=QMTZADthQ*%L+b-^}Sbg+7>aS?ijVJ7iKWS#E zYz-@?1pNxae2wV`nrE3SJ#2gxLfz6&JCT)wxc{MB&L+?>X@x z3Z%%&U2#>@x{m{S6Zl>9JQL(M|C-`!z5U?A0?hyHo_mc)`G13byK>Hx&eZ5UKS|9P zBEOecglxmr0{ueOlc@9--GGC)S_i#R+!nl?z#loEW11!{-FPPJN2U$%lcOm!ZIKrY zK6Y5tv)zIVK|l8rAVruD3avd1gDFZ-&ccQVeNxglTL42y!$ZdN8yF#Ozwm51VppO5 zKE1C=!EL8(^#+n|xhKNN@e0w{^VW|GbIQB>2>#DiE)4m%mEZba5yBscQHU@hbZ z*sFHas2_Ic3RP?KT8x@W%`nZR8CXhc)F+*LfS|o(CLp;H&$dz1b3WasggE+mzi_Da z-LXr)Qx>$uYJOH)v}+XYnK?}*q7%%_{VkPXxB%@Vl~$tU?yM#?VlUeOZfklxerelr z5bS#px0Vr2=B+?!(8&{{{mI<3qGp31xw?TG_DQjS|s#RBhd@D+dx`x`Ovcx!rw47^1*Z=x>zu z-yD)}&^rc-?`7I4EklW47vMTq7zr1V!t@^Td=E)vY#gXdhTHQlemJB^7feDOlJBkh z5IjtcF03SlKsWNUZ6arh_5kq+8b zS#s+%wHM6vZVDGIIZzDKRkPWIWlfJaBgKyllY4bq50IgnrC~zfHX$x6FAWh_kNudL z5|1aTt5E`3vDP{wSV{2s))UAbJnaHEKBR_~@NYlw>PN?M{&|b_7f)6!J_`7ztO^_l zocd!)e@m;8ZDg?4M2EGo0Y1`M)Vjnyv^pohpFbS%^*ig^@E3k+&qA1fya?jAy2$z) zR$+5JOcf~Ud?B7^rnNp8c=k_WFLrl}lz_Y=MI!dI_#*L-Br_GRL$c6(8oSbKWW+3G=1zqdB>m$Cl_)?Qi9T9zBUyIx03`iF|7sCEmG-}i$@ zI!k}D>bFi=t}Jb6xMLa=h)Oi%eD<%ngN|=N`pQJrZQyG_S{#cT`}IV_aCu%!jGYB? z7;D@hzWl&l&niAFbx{g7_x>!UJR*t*at-y1P38l1M#LIl(_&-{a?&X#E6+_XXc4nl z8QQF6TzFbSBpbzHPq!_#Soy;xQq%gy1Bnq=@;gSarsSCExf&Jk>sjg{BIse;&vDSb zGv8NJ4tIVc;RH#VQ}^dcF-*5rNIcqN7-eQtmMr`m;xs^#f2(6qaVkwdWq7vvvK4 z?GCt?a(~3Z)R-;+nCM+Hfun_j8#WBDeetifh|}$dn)3&6QqA;-KzFIg96yEmj++KA zj@;-9;-G_;8dcnxn#W9ty@k+j-!N+IyrhToK}}|F0;g@7=3@sX#KF8?ElXlkxPo4s+YWA-;689qq+J7W z1RW$tABs>d#5s}?(xIwY<4CPN#WCj~t7lMW3%$#)IwbB%XBi}IfpFko%p?~xfN`e!8$o}6MwW$L zy14ZeX7ZR0<8xI>0RcQrZ4B*9G)z4;6Eh|;5dnN>D#Q0IwKwdWKE^=ZM%fICyJE*dg@dDE&5u0KZ=cRx*c) z87iKpc$XQ^o?$2BH*&+#LnjZ(fCXAitcO7W+^!`Dzpzm2=26 zyF-gKxBu_k?sA^hBVhJoij3QoAc|Uz3FjVObst*0vqFAf7h4irpBVmcmB5j(t+30q z&o0MjPbU8@&y~7|TVlN`=ck;{f(Jxy;4*Jagjp&-fQPMnB(BDuv%ZcrDU8)&B<0%G znf)6b4U?t{hu%~Fyo1D@qc&~QuK$aLLz8jpF$^=p@-ZNfcJ{lIk;6I2d>tE%{2X8U z2q6939OUjHKVoH=YAndP3oYKmO@IO(Sh?pgm`CLr7Z7SRfHw*$qCA5#PJ< zk8H>+kUNNf&~f$?T!O&M9I<-aM39k)H~m^_HujOVAV`6Y4V`lPyh0gb$Q1J|=9ZCr z?HB{K01N#8{b`o`qnAO7I)N!Sau(G%oxZ;fEc=ERJN@~g!B5Pk+8b$HXCB_Vt|s#P zzY}`j>oe58EQf|0$~yYK_#kCcZn>`duJ~&zD8NZXW6eL1tj^n@49pD7D3E9T`!(16 zljFrAUQ9T+t*`Iz+H>1}wX(e|c`>6dD4P7b<9mBueV527+p2F>rvohVENz95+ndia zQ?|`sKT1hm8@M~Cgw1($1-9&OG0Z2a*y3Z_aQ|m?yixQMQAXXLZxQ>{T=8q!u8HwR zV-_!h7c~FwSZ(3%Zt57;*8Y5eDSaXnplIBi(kHQ2rz$z%Gd|5;ZHvMy1Y@^$ur`dO4X6Er`s`2i|VrbFmwv@6X3Widq7k)SLCM;dae3e3L)U(F#@DHPrFK zZ+pvsdxtZ==!@g)_d~DNb5bY9XVtC|x-sYDAhNTMcJz1A51VPXGGaVbEFWKAlbh9K z-+XBqWc53wvyE?L1CsSK_nxk;C>Wbf^!YeV< z=EHSeFppBiJT8Y8@FV+4?pmzC_j8zO&>Ilb=X4WygIzfHf%oJ9=Mh**I#>CNu8$1f z517dT^_Ti9jlxDMumS{WpFfQ$3NNU^8HuMX6Dq~nB$h{q@7#~|=R543=8(z+`zP;_ zwK#nb%HT(5D=oymA!x}MQB#v%i6rH?3M}s zuAGy=jO_ev=-?{bpoK(kB2zQcx`mcw^_pXEU~apM1#^8_0^9xm_-NmUzAVuWK{O|1 znkr~@vz!;OTWOoO4`GM;SIkh`d@j9(@<$o$#^k^H$%mqv3~9jw*BYXFP9Fhufr^+q zfqi3`64mr>Pp&UFy|J&$7RS!)Jhcd7X9*CHA;6Yj9XShJ_s2KW;khwDw<|1(jS0jI znjCB6mT_-9*{NB5k=5ZPBTANvmi^Dt#Fy%WN55tY8B)!pIN)e>I+x)1z}5{$4*3w$ zTyH_TM07OK9{@(snLr-GKFS_g9449#LiV>QPoVPO7uUwj>&xp2VfP^cvAp!9e6mg! zKk(JJ-LK+=Sah%ZZ0<0EuI03!&&DW1~w_h0MDVc*26ZiroI6tEZVC-)R%| zOV$6L-bA${qfF7p!0D~qA7G0a$K$XcWg7r&z)8m72`I5eRfjuqe^cMIx!C7us|xDJ zsy-0CExN8 zZmY1TE;UF<{vZMz$=uh}UxSC2k_XFQj&?`TgpZ ze>>#OC1}do&)8Dink5lr|NWy?%0%FkVRx>}3H_Y1?{)zcIkyR%_@Yp!2-9FcGAt{o zv1&OX9DNgRhmdeX@BI=au263`f%1V6_JH{}b0`P1;TLtv!V0W@jy#D)9M;4uU znh{N6T;4>jdLP%tYN}mlDc~wQ(JL+u7b6%hT0}*B_tM=N^x z8Qyc3<1=h0^y3!NQ!oN>dVOy%PY?K}D*YX=Tz!T%bQsK`13oVIl%bZlq|sfcldwTD z+@2ogscwGaogD2O@KK~|z;g`}2kzQT|Ac&Rs$he4caz@W7k{cJoH-huxr_d*x7_C+ z`2o91`E5;o(#p1W1s24ORZp6VeZV@c`E5^s^(N6spyUdDKTG)$d_>OWJ~OkQ{5@EH zLpBLZRt+R%PjfP4s-G-Pdmrx1*S$Kacbw(?gjJ1f*i^JZN$;903z^4p zETZg|bQmTvpa@DKlp#n8Vzqun`6mBF45{-+YS_Qyi6%H!o!XJRt&J#P#$9@_?C&ab zXkuY5%KiG#mB7p5_7|}(Q7y_c^X|vSEumC;?P{E1#40#O;0sZD><@33;#@@ey*Hyo?15I@*CMur@uP4gSmwCJ@- zl&?~?(w6R{BL`h9>k-ZF2O%M-)3%Q_nRo5*5AI!%e>^$gO_vpHfhTlAt%T3A+$v6D zmnE7SIqV&jUlQ|3`PT^13~kk4&HcU9EbFu=@9pfe2Ea!2V|?a2Q#;LTWqfYeaL+bnWua@9sM^HZuy9x?9h$= zTc|FunB#*;dq_Z6Una59ORsL}Y1MsTUqCvGZYb2eJ+)(ct}KsAsJ=ABU&cD+cbCll zYmujf+htmp603?aROoLhKA>tB3;t6Czr!0U9p*n~N?xupns;p1dWKA;v86I!&xh{3 zIGXw1F@v<5xhMTh-@E=~*+hu{)$*7QuKn^tKCs6noqmy3sNw!$$B%B~$#2p{-75>^ zHuzjsO>T2*>7wl_{RUAAYX<8CSG1SNyLXFR77peI|cd5ALuqi@6(KMs_qT6$>MMW0t zN(klGss-J)As^X6x6%m^6Bp7=f46faM)23!&9IMTg{;=aE<#^Mb%F#FLxfZ_ZhWLl zw;(q1Tswgg`8r3W3HWsV7*~+4u$jY(0bX2~_YO)SkE4Ih%4k2f5-rnMp9tPKRy0um zs*GhVHJ?|ziu(xwA1<_w`>V%;3Il8had)Rd$P!Lkybk;+D3uvB`$c!=4Ufb;XqBFX zjOK2t?|IVc(9ZFeV&sOjHVaRjc?;?ThGXv|D$4Jg67zPo4R^#grg@<%0-Jri&(`@= zzU!bf24nAveOh#8w$xi*=;))C3)Tv-TICJ|&S@3)bi^~gA@>os^NI3kn3o7EhZO4) zMoru89Qga|nmKe-GhA_XvQPH0;|8OFxqyo%y|I+uMYd>_0?-MGBY+0*?S15cW^p({ zrJXtU(n48zsj3mJFvqEL4DP4JvwDoD&KeKlU!`t_GiV0T1ynA;K_fZixi+MZ;y4^> zK&u2+l6F~6{nr)73g(--x&y-76&Yu{6OpPuJj9`(I*YN~{(TcyAgn~fa`*kjz zOk^VP8&dGWm9(ch#38?e(x|+;n0>DdOPnsBl zAbKXG{)unT2$*x?3^s-BODolg)q#l$#C-rf!$GS6+n%Kr14Lmxfp+e0-@1m9;%x=; zH>8I+^6yR#0ExKVQ?x=oWb>8D*c~x26%;!UcIVyF@$!5CU(AnR27Y>dRNk^kk|muC zc;}xAP<;|u-~K!n-+ji)w^tm~EZ%+_wMvVK^$2QkRQA)e(L%dkjwih0+%sYRlfV;F z38t!q(~LVrPWJjlu z$RFRE`5zUtfxt;$?d{)3EB&zZf4grZD~`ZFcjbH*9j8`%d-iL$%oE9i7tOC^CPnNY z*ZhrsqjAl?a>jgW0&yOUZA+kwCzCeoLCiy@U_S`A=9w`-_RSbr+lR{%2G&l0zW(I# zcUexESZGC>)|cI#Q$GB`)5S92VPDC_il{$?tS(x;MuSt0q?K%aaxvBrS{ykJNVuAa zK3?A_uIHrH*ZUpb-NLRzn}RbmF0Dzjz;(6(16zV7wX2#ag)y~soUKZniMBPQy}`0QH`ms!E7&GLhV zlWeY_6+=7|cv^Ohu`)d{sxR_;Fl5jWCv+aO|32l-<79 zJWWtWM1MxKSe>>~_i(zgoEAT@4MUAh-BM%FpDYrYtCsbsH z^;as~9ZYIJD5{f-Nn5uA3b_`@_xsY$jHz2*bumb>R#k?B9o1_Nx-L zx!u=j`SrC?8*;FF%Tg|%S);3*t(b2}ZdqKb;^3(_{1Monlb0qUYt=qpWI2R$C*Y(f z(8;`)TCH`#caCl54CIQ&yjvOPv|2y7vR8<0ZA6YPc0XEXoa{NVMc?K#?JZ9UXXS8B z$0;ulF7{6#Gx>J8OI^(Ue-L&fxDOnNNcj>$T7u^Ni_lv@7*Du_Q{zbrI5STQjXO2#edr11J<1*E6SUcFV|Dd2f5vy`Z^*Qdo7+vlI}a441xVB4aSizY0-rl zsC8{Wh(@xXNPLOqMo2aIrmp72^(Lv{n*OD{@#EixR9@%zd$u)~eoA$Uu`8^K8$M5Y zaG{R}6-5nnu?Y2cbW(10?AvL-QHFR5qEsC$^HjB9r1+Yl4zOj9b#RLtxAgL@74&}D zC{5jBT<=6f*Ov?e?+OQ%chgt+mJ4LDQ2IfOYe*mOJZMqZQBoDG5Va&A#eEh2%l`sg zZ0(gO%aLT|2Hq(2L@YDP(o`UGS2YblNoSm$pp*L7QT26ZBjos(M!Dqut@tS*)L(fK zR`c#8BP=PM4I@95;)L@)GZ#DBz+KzHErxlxln_|eT31Q5rH*knqf~5c)ee!Nxa9>0 zNBEUtA&QA4Y?*@YP8)oH`*r-4w|}bKfo&(Y*+0JVl)EozK5lK#G)C%(@;}sNIK1rm zAA9GOu|n#&TkOCrPI)^EDkH&jC{s|Ws0EIc^80Zd6-pVVl^o%Mz+G0`xr?zNEvexw zeVu@oRcf0}dK+(sL|Pwv0x3vB6!;R;rcs>oLyTfIqkNV%X%3ZiJU#-p$Mo%bZph1R z5qy)EIj9*um)NPy8(Cqz{bGnn?s^}~sMIJm9m%OjOOO_OYvi9oebyy*^j&=31$Srv zo(rG}=a5$H_sKK!5|tr@mgd>t0n7$AR_}T8xqs{KBVT-%83zwajH zFA^pV!Vo_Uj~r5Y_l^k0A0D{CjB}(wjaEUH*xkZsgNYB+hOh@Ke)BSKD>aE}@|e4Y z{lC$dVtWMjlP;YVrBUrPV+)m@HS*{&x_^kxKsQB&$|T4V73BY(KPZldrrIY|hBTCKs>Q*I?Shz8+EUM^DoE!r!hVOUa)st@YT@ zYMU+Y&dsSN7KsuWV7(`Pi1wOn+}%!(Y^++$t2)b!)j1wFu>01c6I=IvqUyrsVE_2! zByBCs1GUGs`xgFNR) zTgISZ2X_#|&6y{xf_^WZ)9TC`%|L@4E(0mdWQI*0^$%L6ouc0vTk9^4%`&~LaF%(> zWvK42mj6aCi8PPP_N<*KWwlZ3qyx#SuAb3KIgx1B|7McQorXAkWWlZ#q0FvqbrEdA zTG}H-3IXyc@o|5=M6@_dCi=>8|ChxMzc#aZpb<511%ZmApR{agx{LYyrL7yGz=&gB zF6W(eG`DAYiJZSXYUz*5+)+c17nlKYd=*7|-hWzzC2q&J?^7?=DzxJA6%9Oa7=DUuBs_Yr#A3z*LyHLaD; zvcX+j$;0H?rzy9=`b1W4ADQbkFiz7J?h}t$>KONu6IsjCSlOi}zMnU>MLxgwH<*mu zIb<)sv#5~2DZe8(=k!VG72oPm~Un6+1p%|HoVj2>A;n{gmNA4ltkw zAES@>^MfV2rfX!t#CV4z6E+MUKT7Gaj4}DdHF~^(XS{7`Q~%p=$i~C@Cz99s$v63p zav5CK^$4_MtpD9r!kE6rW|YcLKBnGMUCybh0Yr3ThdUT*11iv=NOxA;xv+0u{%}*h zXMHg8JvMIXzaXaMeH+E)x?Hy(Hx-}KBPy3k{`JVO|Igf^gh1PO?6MNGFwe)v zYPEM@9dE4HhY~5+k9mHl@p)FCbY85aD-_>6NaoF|LJChojWu_)L0cIo!OCM{f5>{R zyhARhXpk=5fUky9|!*4*hGz{ z`w11qw$n_C{EsnaX#o#r0^$&;8ig>j6$;a~jTsaD+6mo`VucEttQG*VD|~}2U(!oT zgPkq!S=Dotr{6FCh-HG-h-zXY3G*CoZf6K6pPaKzj9Jl8){o2Q8-&h;rgLhV2pDaL zwWaF@NKMHuR+p3Q&2igPF_<5A0>{8Z8SEy2%FuCYbTp4@%;mgCxgTrR`46Il1D}K| zl$nX|g3ay*0p(S7M8=w_9bQiGu$|sDu8KhpKEz9prXoSp?&%%Ts4Lsx7<>@n{FHZ> zpu`noAcdtm2mal+fK%?43h45S%n}4!UF+VI5Jfmvoo8j4P1LJkv{v1WCYqLJJa9NK zL7lmacRzP*?kCVy+yV0t{#Y1zWL|J#oRJlrp>&p1;fk7<+Yrf*)#aZfn{dOzfNR?* zKEK`Cc_00G9b4d~V`7xgjZoUp;uQU)VMh7>MPlH`>NBJuVj0^mXuS|N2;;%2_b|sl z5k|9bW2)Z0H1!2#Z|Lj0BmP;`%gqk)QhFw!@uJoXo@<>__1yPHd~XcJCSBdJN)rg8 z&3(iEJ${ZT=QjV}{C{p((fru2EiW5UH1hEPJAJL2R110lH%-(aS9LLWmoB70;I_!Nx1rO13k@_2hod5QyG4<|8K45ioiLayeo z*26Q0x#(XH9!b!vy1A6shr52(%!@h*Sb>%%+XM?Gy^g~X&vM1BXid))I%q4uhEzIg z+>b2~9{qQzecDSRXdqnes1@uj(bz+p-++*|Y1Wymn*6cB^x)POPDWN3Uh9%c<3z`B z%&xbe-S5lj3YJxuytaIIn5_x*>jL$!&biB@_zt*0m}wgi6#*4+im4duS=Gj(wG1kL zN%1Hc-63zYjuQ zv;>=zkxLm8t6bi(!HfWiJr5N{Iro`i{7nv@^8U0H?VZc%_BmW{2R-JR5q?cbqnG3q zmT}Yq5zE6|A9bC)Y30Oo_n5o1+;v)f!}e=#RWW-~ecs1!Z4m-Zlbe&jC23%jv_a zbC&is)-~H#>|#iM!bHO*groO;P@Ti|<#MkeP$y4DDwX1wTh2rI*^hL@=a`8@Z_WG$eEwR(~I z>g~>*EiCdqr2Z5Trz~6&$`bGrvd2U_oLXj;C;DoC1EIhd2{40*fCD<~ia;s9g*JvG z_A(FkSJgt*mntx0TIEg{UVj@YBJgqTSCvsM%bwY$3vhCE2n<7G>+@=xhcUxTLc)Wv zPVj4xjKhU!P#rGybO=pEM8h^k)UnUMibP+yalkRXwPh$WNNeZ9&+4hWk2Mf%zFE9%QAkqHje*Vc)-R{_n)JSjgEZAc+^r< zg?=axON^g_wiaZsD^NXt32Iq#Vkd4R*`~qVJWlt3?xwPb!XR(NW)IZYUG$^-^p^i3 zGvr}};qOF{JO?7*cXZ9|WuyMoKP{rbi=KZ@iO__Xd4a6{@dZY@F2m$1rOECG_LUxk zXE*ZWI?rC+Us2nW2Oj{D2psJM`j&^a=GY09cqU7BDE8h?;lezU$g`|%-_Va=Xvh$eSdD;{Vom9p~@RI@B`;9;Ec>6J@s*M zKD*(uOLI~EDbIEszdJReu0V%Hcg6kAY~nj-CGMbxul6cyFi-pCv3@z%D>$ZI(x2CG zK*Y`C`(a$V5oMm)E~U(GkM5#HvnGG0?FzLT_HgUR6HOH9t^@ue=I;&99}G?0EVzGw z4Hh|gbh1Duh4y;z(j882a}&4Vm%lq2r3vsy{QTf1MX28#JGCl6vzKA{LDqsEpeP+Z zba2io_P3f25=gDw(#o??XNILmR$xYV%mXoa>rmQBG@{13NY@pi}we=NN zS2OtZU*Lwe&?~e$!W-|W+Ge`~-2$!?d#VvD)b)ua2!zE;8 z5gR6(zLcdG-^e*ew1wA*Qq98A1xq1Oh>@!+PbfNq+MXxvBSdNNt#i&M2wk~0HzVoV zBG7>Xno$17Gbgap@<~o9^PY{pGeo^MngI*FwwROWV!j^{h`2it-eP^-l>K}a(GIiN zx#<}bl5CL#sgUtwJ$G8#Z(0K_6?<<;u%CAe%5L&2$U`q>&OUU853TORa7S`a;UGQK zLt;R=#RNa)` z7xBrO!@E!-`~5)&4dH_1gtUoHO02F&)XNH;W!Nje*>6P|z%5Uv#>qH+ML2?tA}v3l zeYuJOQk6VlI;rs#;Cmr|TKnId!A~*?0xzcZ(z3V()6*e`eJaBy6QrWkN3)UVk$*YV zoQX`HHa90FjK;MUsQw((L^qK+z1O`4{kv2t?Eb0j9?wUDpnjW6nE_T~eYzK9!m7G) z)fiU3tPw5U#m@IPSzy;Ro!C$GfXMHx$(hkk&jfRgAH%6T8U2k7|AmIP&RX?6fVI zwYB2AY7Cw12~S%1g-mFgI$B7wdVy_Lc47kwU$R90J)A#;JVjHh(k!a}>U#Cipy*O9 zd{c2T_ZEJ&8)ImDd;k@U{hGlD=cSgfq3CnRN7cjogA$UxIQ!3AMRL0-Ku@ed@X?E( zRzW`&S@J^W_%Qj-e0ue7F}xAYNv4NmfOstA6W#bCuw?>3Z^Gt|0`?k{TR^y-G&Y(< z&pqJZ`hc52OdzeGy019OQxzf*!PTl8>@Niwz&dXqowg|%q!ATTC(SN#pz|JBZ>lJ= zgjmLWvn46Sj#iQFlF2}QCI*y(qr*xc3l9-uFG)>k(Bde@vDBhRv7hbwfWU+c`=xfR z0vrvvZQoK0#AOPK6$e!RdwFxfT`9lH;#y7&AP1t27mE)C2%r-5tFatbtj`~B0>f@Bj zr+b`rxH2T}1RJ29-mb@m!)|LIuI@d26@!}fGrB&7n8F9H!8W;6%&~@y1BBv>=*VfS zW{+LjE}$aKg@vT_$j5%m^|4uG|7aYx}9O9R^FsDt97`|;DP9KbWLMaBhUzibQmi!EPi1$aG z)Ou7Q|4DDGV&)lW1tq0Djh{W&+@Gi|!p`#(MfwyhO6X~;z1=<;h zW$lkUVDF~$Cb-&Mnx%VA_NUx8<-ch24T>Xfs zXhc^K#dV4~9#>VUAtM`OEra44O2I3~n8Va{Q9^(K7T_&2E}upTgCdKleL=G5e$H_` zw&Mu1I>%%!eRNdl?(RK*#6EsRM<(_gU)$4kB@2f1-{`q)Z}(6Hsn5e1pcSdv6Sps_ ziRfgFR`r|*T}roUC){DmpNS}H%N0;uDm1NUB_PcFJs%te|;*Z`I^c zXt;$VSvtELP|VpKreNV&VUGFl;|O_o^IZa2QIWH#ym*$|Ozc|OS6=#WuTR=NsCtqo z`^x2I%2H)aF98UHz+e0(`n4oy$qjS(|DRU!A119F^{-gGf4hEEv?%(i*X*dy$?Q^e=;tSCN^0^7b(o8_6PEI7$Rzw*^gBj3gq}`Q^x&$bV<4-1qCgn5wN#t*eydO8O}3H3-{W zBTXz?!hdfsqd{uD*k$U`Ze6I6*r2H>Not)g(N5*dC#+YR|5G9l<;(f{xSB+K z$ol>xVGH6oDo83r-#^6hMqW{*2UT1CK#UK+U-`mEHCJF@Z)x%G8O^?n`eF97LEcf> zy11RPIEX>ktY6DCm$Vr_O@^13;3J(5FH}>OX}JDLVLsgtA&L~2pU3~eB$pnNW7u3K zyDdZ4Wq|W5B7e(x4~tO~riVF6p3urqyLy39tCj&l)u;(CO%D&H1?p+Ko*jg!-W1#< zgraaJC{D`5Fy7TiO;<0hYRK5XqA4KTg~%Ex#l3P^?rO}uY^PB8zzN8YIv| zpjlYF8}*PH!3qAj6kcb$8Q2LyP~?n%y;*-_I0%}4yr=uq_awDeON3)dnk8*gHH>w) zS0%U&T53fxv>>+9``|mnjQ|X2lHggc3{u;51Y08L!!ATbuK;)!qt?HdNs?LF> zj%g%Lc5oQJGX=CiS1HQ6Kcljo$h}oses)r%K(8dnMCx$(7C--}!yf5`c<$@M07O;v zvu$PWc_FjLc0~uuI=R!+c*gVq>|e+WihRI2n^je(r}laaIVwi9r8Pz>k{g-=-^x2e z#!ToSbRG|FWlJ#lQ}SEu>j+0_Wvn-O>{<0FLOxZfEGxZ=;8p^WgKSB$4Io zC^3LMd!Qbk6U8xXRa|0fr|eXp)k>p3emiXYoYg7xd4XB%lEAt0*8L61+Q}f&WhBeB z@Yn25!Az{uDOy$tbW)r_-IisSH%uinE4LHwIN<`S-L-~|&Y5-5e^>fe&*t1CD}B|$ zrMGyN%&$a*TUggPDJM%EV~=Xj7byunxif*zVa7^DA-yEm*c;R*c<}-CYbQu$A57TS z+Fk$&?tdkA-BXcriV}?m-m@o&sp^ADw2Gaftv^~;Y1O8YV=>`Nupq8qKki!95L3Xc*1;P7pU}D3`}+P99}~K_)*3drS2?C z^%o?o8}HZ^++KS&al1N#d~$n5l-(l09XrVIYGVcq<){>n&UHNiAtAD8dd7gje5fl_ zjC}L53;tLsQH_8^{D{brM{Er2oJQL0pNoZUL(qRBR18fb#gagaU+Y6pED8kZF0XQb zY(ac7YH;m$X_TE<3(nQQ#OteJ3q5=Air=P@7CGDv2i#IaSa|h?m=?HAGHWmGPO-go zR)uJ$Q!4LMajHZZ8CoHEhjVZ6XqWlth(i-edrKGK{;M8~km>gy@$fhgC7NENG}}d! z6_O-hpf%P>=)ybWcna0?*1#gGVx!y8GWX#oR`*j`(&P?swGbUS~V(W@^F* zicC5=4|x6TZ9_%q&vFB}C`v;ksSYIH+$p=I%WR5?hzR4grCTzmA z=k$2Z@{zNDhLJ;a1O%9fZSI1bd;7L%*zq zK_4)jtz(gY8XO2NX}?1L>j1ks#=r>US~#>7_9(Ifrb@skQeFN)yX&VCloxmNe7Cp; z?5*6mFZtQ79g|o)RKnx3C~<13KBkhG*$RG5Eg$!_<=AiuLakZEE9456H62Yo`v!Z0 zbD!xBjyQkv;?%KGy%yzH@h9deXL=>3^LE`Mvy$m$gNU~`!3l`(z{TeF?lsUvwOk~CJwc1PDSBL_>;f7$%Oe&06z zl!=y1Pjih7R`$cFo%G@Kwm|_uV=Ts;&e>~~PlW$Eh}P+j6vv(LW}dNkd=7(G!xy3G z()|a`TUUBpDPY{W&2V|(*)aR&X_4ShEK9uEnmV#C)E0ikZ!#Qoch)XZ8TPPv=La3T zv?mUv|L13GKAT4Eo7fUi+B0sx|JuRcrU>o44f3=EDCLv7>H9T1TI3;%VysvU28~&? z?`w4H0Qcebdn>@~#0WR-V!I4d>ic$z^V4jyA>*_cc8E~CbNc;3SSavj(&9+%rZ8tR zKB`*R3*l;_yOv$R3Z5-F(3o{X^nFs6RcJ?fB}wihwPZ%V>oLUFjdK`Fk^#?KMQdf) z&`OBR1VXr$ekfhi4`cuuXUFaQea;H(%ywzJ+a*ei)_Lx@`%U`5_2>OwIcB`utEwhf zoDA)+mfmyNcYbE>oSh<(G^qNM?Ok)2`XEU8HuFZ6;AeJj9$(m>ub(qlPQOzq@M2Vgy5mMqU0>!0T15S&|`f0hINeT{xycUT&QAB37p1G<9lCV(ULGmHqPI z@2TYGWTQ2dM99Gildc(gupXW`Q*Bh_Q@;rL!M4sfooHqXG8Ng6?-Ht=m6)G*o@Gs2 zM*5H1PS^E^KT8{@4%{m3rP9OIx6ofc2g-X1F>SVv@jtZMK0K;-GVl=MAI#~S*4MtU zJV@JHzszOx9KYd3h>*-w!!Jo}%u&u67ZaAtpY%KwXfgili@Siob`TGmJqF;~{7|ri zX^~2Hr*_@51vhLdZe#sX8#m5`4^4F++C>ItjaF~XI^#6hBYM;yb(|9DdV1)f@y*w> zI-jGvQR-7_(#y4Qz=FcF?*iLaX52?mzi%R6wN?mhvvDQpuq5QDjA-r_DIRIZvJME@ z_N+Qpk&9!bs)N2jp~my4!ZZ|4`cENOO~kPe<7XrLPOq6ZPcA$#pxo~Wn{Ma9U58ED zCULrN7Cx^6YhUPmSjjHw0mi{5=8!~&z8!~cA%uv|4nBnQ?BTVdwS}3)y470XM4ZN?9~d!y(X!J3B$IK^>kHjpic1 z!YYdI|5DeRSk_e&GB6)!&n1b8q8M8JdhgGmY{<8rY?6#+L7qza+>i%=!s^-P{boHe zmjkn^Kp-X5B{zj2`tZ6%*Hk0}_NXE&+xJM5lb%CWe{;~) z#IdrKEEN^fE#_^=wjwuT#(oP9*x3}PTi;9Me2sEKY@I_A?8olgF{8xt1F1 zX1988b?vbn$VaoWpkD>ikh5TL;*cB>0Y9V=19`Lx5b|{ z;lm|w%d|bmC2I3??^ZQA+~*G(%l5oJ;G&cwQF+rP*u@b`S)X9R(%KNrwrT=}#L$6* z(&LvjJx62ooQZ2Qs`$6ypSJ1{Y-6x7 zq1MQ>&gqbOEUMU0|@*eCF=8#?acHySOq0YE-J~{ai?Z5^dg zB5b1>5D&{JyD%ew5oKo_KyO3&yoxMvT&Z=uM!z7!cCCE#k-jP;_7M+?Zpybcj(!bn8(%MNoHOA+iz7w;K-RX~Pq0qA=!OuCa68 zzm`~d{0X{Vfb!D$i5WFHebA&{K`&I-|K#8wA!D`3F^fFY-Irdu?RqK>6W1H}zC~V% z;?UA~@zYZ4gTZYvXfZj2iPyPiamOoeH4(COyX#SC>NBw|Yw2N$Jw}`b8^n+M$In(5 zV*&xP2iALTlZ#A%BsE@CEea#@M zAsJ1OuArG?63cN*pG0X?UR5r^)Wx9#B@6m3b#?VD>}kUl$mH!L7$_MxvqFpZSq%jB zhW+?}7oWOxlqfsXC^1>Z#@*Xe{`4CC+Ud6s{HRtSp(S3_S!Os!V+nO(gs1x2M@arm`q@zN+Q+n^H@g9<|L&Y4Xjib#c-yG>g ztVvS=5cNJlq*bx~4P(^%(gh_M1NNiV>2FV)@wc|+10qWo9?FsH8MaH-7mbA3D$R1p zFYePSPdQ7R{hDbV?&^)VYKGaDv;al2nRayHS$4Mg@7&=L4;WSM?0!`KyVMZN z=mHP6g*zbCdb!j78al=>apM6yV?0V*kY=K}pks6xgysq0;L+qGpmnr9 zHMbF0pNXM3RxNYhjsaW%?{+(-kp0Y)teCSAu9vr(2``;NbwO7bAyj&t4ajPIg#8(Q z5&yVfJ7t0lHQh&_6XKjMA$Zdu=Nxg1c_Eq+LCeDl&|A7AVFKAGKUEVE* zp`SI_ln;?T9Ry|%2>K@n-PVbM6h_01{@OPRb1GgM?`7@nNKx|xQ7>-J(8{J_QT~e; znGsVrt-AnBBGCV945=<50G{2^J%y6K6}6x| z|KAO^Qm6;Z4q2rBVl|%d#EcpuXv9(8DZh=qsu3#8-uR)3*hOn5GOylUJ;gZMKdW#4 zAi&`m`(<)j-f!BMGa2Aee}yBQk^QpIfB1J|0bk_RJlSVq)qc$VXsD))h1LjSE}yaL4!K2;xK0vcSL)A(VBJR1WTa;vE}qD z&t!=)Zn-_%EV-9IeCtx$F7QDd)(`k{2>FFFJ;GjFKZ%F>&Dr4jWkXfAvEC_PAXYH& zsC&~3`Epc}``b(GA#7s7Udo*|(_opr$iJxpkmKcOW*e-x!Qd5OMeZZ>!UC{#Uf{Kp z!O3~sfA#0tzk5807a}ef0e2pju>$@ntF#4|IkSgX#r;Oorz0nxAWDDgJ42SEtGfwZ z^r-KKw!mS;P430|C0pwIXkvRr@a;K#cR{^d^G1O=-SFA2?WpU#zZ@O7oPQ*b!>i%B zPxskTTqCaeeA9A(TYdeTg>u+L-AgiMpSoHwwiiRXO^t7B3s48{?hF?mOk&?56|uaC z@-IjBjW}@sKH-1JzA$?(n-%tM60#kjddweHC6uw0WjflnkO;WeXhH2eP6cjzvQ)I< zLals4k9d)dG0ex}NxxKF?RtKD@geL}?=sqF;+Y6M3KYzUXZx^FqUq*R7-XG|0r79$qi( zDOeqCzcwi5tTeIz_&JoVqbjGjBp?$Rca@_+4x12C1uMTS_{qMl0Xj~y|4n#ch6P;y zVVoHf+B23#M7i@t1FL8!iD-T)Lz$sfMAeEtHTCSZqkg@uOSc#8Q@MH8PlfdzA1*us zThKjAETIZLH%8rq9=^|YLX>im4ic;fW~m$TRI3*+8g2?engoM!%2{Axh7R4ZqB+g@ z8n$=Nr|QD47s6T>lVdKvi0Z=}E9$2MdMA`~2LW^9URZFaowLvB822H1EAb3HNT z^Tka^+Poy@?duzP*wI76rZ|Y^!j(sEvO-4%Rkji$Cm3EqJxjuL4;Dqkr(1QhdYyqA zqc5-c&Z7(}1eQ1-Va6jc6qjbXL9E7TB&Danp*h^2EE%he%6Yf+r@^XZU|}yK2B*{z zYW7r4apTawEpGf?=~l8;cs z@^-LIBtsQHG6q=hba=^LNKJ*o4zEK-Q{?avYxauFm10_hZ|ykhf45#ODC)W>!(r_e VvjYGvX|+D+= -1.5 && pt[1] < ncol+1.5 && pt[2] >= -1.5 && pt[2] <= nrow + 1.5 end + + +function draw!(g::LineSegment3D, 𝐨::AbstractArray, 𝐩::AbstractArray, col::Symbol, p::RecipesBase.AbstractPlot{<:RecipesBase.AbstractBackend}) + x = [𝐨; 𝐩][:,1] + y = [𝐨; 𝐩][:,2] + z = [𝐨; 𝐩][:,3] + Plots.path3d!(x,y,z, w = 2,grid = false, box = :none, legend = false, linecolor = col) +end + +function draw!(g::LineSegment3D, 𝐨::AbstractVector, 𝐩::AbstractVector, col::Symbol, p::RecipesBase.AbstractPlot{<:RecipesBase.AbstractBackend}) + draw!(LineSegment3D(), 𝐨', 𝐩', col, p) +end + +function draw!(g::PlaneSegment3D, 𝐩₁::AbstractArray, 𝐩₂::AbstractArray, 𝐩₃::AbstractArray, 𝐩₄::AbstractArray, col::Symbol, p::RecipesBase.AbstractPlot{<:RecipesBase.AbstractBackend}) + draw!(LineSegment3D(), 𝐩₁, 𝐩₂, col, p) + draw!(LineSegment3D(), 𝐩₂, 𝐩₃, col, p) + draw!(LineSegment3D(), 𝐩₃, 𝐩₄, col, p) + draw!(LineSegment3D(), 𝐩₄, 𝐩₁, col, p) +end + +function draw!(g::WorldCoordinateSystem3D, scale, p::RecipesBase.AbstractPlot{<:RecipesBase.AbstractBackend}) + πžβ‚ = [1, 0, 0] + πžβ‚‚ = [0, 1, 0] + πžβ‚ƒ = [0, 0, 1] + 𝐨 = [0, 0, 0] + + # Draw the world coordinate axes. + draw!(LineSegment3D(), 𝐨, 𝐨 + scale*πžβ‚, :red, p) + draw!(LineSegment3D(), 𝐨, 𝐨 + scale*πžβ‚‚, :green, p) + draw!(LineSegment3D(), 𝐨, 𝐨 + scale*πžβ‚ƒ, :blue, p) +end + +function draw!(g::Camera3D, 𝐊::AbstractArray, 𝐑::AbstractArray, 𝐭::AbstractArray, scale, p::RecipesBase.AbstractPlot{<:RecipesBase.AbstractBackend}) + + # Origin of the world coordinate system. + πžβ‚ = [1, 0, 0] + πžβ‚‚ = [0, 1, 0] + πžβ‚ƒ = [0, 0, 1] + 𝐨 = [0, 0, 0] + + # Initial camera imaging plane. + 𝐩₁ = [-125, 125, -50] + 𝐩₂ = [125, 125, -50] + 𝐩₃ = [125, -125, -50] + 𝐩₄ = [-125, -125, -50] + + # Initial camera center. + 𝐜 = [0.0, 0.0, 0.0] + 𝐜 = 𝐑*𝐜 + 𝐭 + Plots.plot!([𝐜[1]],[𝐜[2]],[𝐜[3]],seriestype = :scatter, ms=1, grid = false, box = :none, legend = false, markercolor=:Red) + + draw!(PlaneSegment3D(), 𝐑*𝐩₁ + 𝐭, 𝐑*𝐩₂ + 𝐭, 𝐑*𝐩₃ + 𝐭, 𝐑*𝐩₄ + 𝐭, :black, p) + + # Connect camera center with corners of plane segment. + draw!(LineSegment3D(), 𝐜, 𝐑*𝐩₁ + 𝐭, :black, p) + draw!(LineSegment3D(), 𝐜, 𝐑*𝐩₂ + 𝐭, :black, p) + draw!(LineSegment3D(), 𝐜, 𝐑*𝐩₃ + 𝐭, :black, p) + draw!(LineSegment3D(), 𝐜, 𝐑*𝐩₄ + 𝐭, :black, p) + + # Draw camera coordinate axes for the first camera. + draw!(LineSegment3D(), 𝐜, (𝐑*scale*πžβ‚ + 𝐭), :red, p) + draw!(LineSegment3D(), 𝐜, (𝐑*scale*πžβ‚‚ + 𝐭), :green, p) + draw!(LineSegment3D(), 𝐜, (𝐑*scale*πžβ‚ƒ + 𝐭), :blue, p) + + +end diff --git a/src/estimate/estimate_twoview.jl b/src/estimate/estimate_twoview.jl index 19603ed..820bc8f 100644 --- a/src/estimate/estimate_twoview.jl +++ b/src/estimate/estimate_twoview.jl @@ -1,5 +1,5 @@ function estimate(entity::FundamentalMatrix, method::DirectLinearTransform, π’Ÿ::Tuple{AbstractArray, Vararg{AbstractArray}}) - β„³, β„³ΚΉ = collect(π’Ÿ) + β„³, β„³ΚΉ = π’Ÿ N = length(β„³) if (N != length(β„³ΚΉ)) throw(ArgumentError("There should be an equal number of points for each view.")) diff --git a/src/noise/ModuleNoise.jl b/src/noise/ModuleNoise.jl new file mode 100644 index 0000000..c6bb26c --- /dev/null +++ b/src/noise/ModuleNoise.jl @@ -0,0 +1,6 @@ +module ModuleNoise +using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleMathAliases, MultipleViewGeometry.ModuleOperators +using StaticArrays +export perturb +include("perturb.jl") +end diff --git a/src/noise/perturb.jl b/src/noise/perturb.jl new file mode 100644 index 0000000..d9aef1e --- /dev/null +++ b/src/noise/perturb.jl @@ -0,0 +1,15 @@ +# Assume homogeneous coordinates +function perturb(noise::GaussianNoise, Οƒ::Real, π’Ÿ::Tuple{AbstractArray, Vararg{AbstractArray}}) + 𝓔 = deepcopy(π’Ÿ) + S = length(𝓔) + for s = 1:S + β„³ = 𝓔[s] + N = length(β„³) + for n = 1:N + 𝐦 = β„³[n] + D = length(𝐦) + 𝐦 .= 𝐦 + push(Οƒ*SVector(randn((D-1,1))...),0.0) + end + end + 𝓔 +end diff --git a/src/triangulation/ModuleTriangulation.jl b/src/triangulation/ModuleTriangulation.jl new file mode 100644 index 0000000..14b0532 --- /dev/null +++ b/src/triangulation/ModuleTriangulation.jl @@ -0,0 +1,9 @@ +module ModuleTriangulation +using MultipleViewGeometry.ModuleMathAliases, MultipleViewGeometry.ModuleDataNormalization +using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleProjection +using MultipleViewGeometry.ModuleConstraints, MultipleViewGeometry.ModuleConstruct +using MultipleViewGeometry.ModuleOperators +using StaticArrays +export triangulate +include("triangulate.jl") +end diff --git a/src/triangulation/triangulate.jl b/src/triangulation/triangulate.jl new file mode 100644 index 0000000..c8dd2ab --- /dev/null +++ b/src/triangulation/triangulate.jl @@ -0,0 +1,56 @@ +function triangulate(method::DirectLinearTransform, 𝐅::AbstractArray, π’Ÿ::Tuple{AbstractArray, Vararg{AbstractArray}}) + β„³, β„³ΚΉ = π’Ÿ + 𝐏₁, 𝐏₂ = construct(ProjectionMatrix(),𝐅) + N = length(β„³) + 𝒴 = Array{Point3DH}(N) + for n = 1:N + 𝐦 = β„³[n] + 𝐦ʹ = β„³ΚΉ[n] + 𝐀 = [ (𝐦[1] * 𝐏₁[3,:] - 𝐏₁[1,:])' ; + (𝐦[2] * 𝐏₁[3,:] - 𝐏₁[2,:])' ; + (𝐦ʹ[1] * 𝐏₂[3,:] - 𝐏₂[1,:])' ; + (𝐦ʹ[2] * 𝐏₂[3,:] - 𝐏₂[2,:])' ] + # 𝐀₁ = vec2antisym(𝐦)*𝐏₁ + # 𝐀₂ = vec2antisym(𝐦ʹ)*𝐏₂ + # @show typeof(𝐀₁) + # @show size(𝐀₁) + # @show size(𝐀₂) + # 𝐀 = vcat(𝐀₁,𝐀₂) + # @show size(𝐀) + # Ξ», f = smallest_eigenpair(Symmetric(Array(𝐀'*𝐀))) + # @show Ξ» + # 𝒴[n] = Point3DH(𝑛(f)) + + U,S,V = svd(Array(𝐀)) + 𝒴[n] = Point3DH(𝑛(V[:,4])) + end + 𝒴 +end + +function triangulate(method::DirectLinearTransform, 𝐏₁::AbstractArray, 𝐏₂::AbstractArray, π’Ÿ::Tuple{AbstractArray, Vararg{AbstractArray}}) + β„³, β„³ΚΉ = π’Ÿ + N = length(β„³) + 𝒴 = Array{Point3DH}(N) + for n = 1:N + 𝐦 = β„³[n] + 𝐦ʹ = β„³ΚΉ[n] + 𝐀 = [ (𝐦[1] * 𝐏₁[3,:] - 𝐏₁[1,:])' ; + (𝐦[2] * 𝐏₁[3,:] - 𝐏₁[2,:])' ; + (𝐦ʹ[1] * 𝐏₂[3,:] - 𝐏₂[1,:])' ; + (𝐦ʹ[2] * 𝐏₂[3,:] - 𝐏₂[2,:])' ] + # 𝐀₁ = vec2antisym(𝐦)*𝐏₁ + # 𝐀₂ = vec2antisym(𝐦ʹ)*𝐏₂ + # @show typeof(𝐀₁) + # @show size(𝐀₁) + # @show size(𝐀₂) + # 𝐀 = vcat(𝐀₁,𝐀₂) + # @show size(𝐀) + # Ξ», f = smallest_eigenpair(Symmetric(Array(𝐀'*𝐀))) + # @show Ξ» + # 𝒴[n] = Point3DH(𝑛(f)) + + U,S,V = svd(Array(𝐀)) + 𝒴[n] = Point3DH(𝑛(V[:,4])) + end + 𝒴 +end diff --git a/src/types/ModuleTypes.jl b/src/types/ModuleTypes.jl index e4ca9d7..d0e191b 100644 --- a/src/types/ModuleTypes.jl +++ b/src/types/ModuleTypes.jl @@ -10,6 +10,6 @@ export CoordinateSystemTransformation, CanonicalToHartley, HartleyToCanonical export CovarianceMatrices export Point2DH, Point3DH export HessianApproximation, CanonicalApproximation, CovarianceEstimationScheme - +export NoiseModel,GaussianNoise include("types.jl") end diff --git a/src/types/types.jl b/src/types/types.jl index 816bca3..163e9a3 100644 --- a/src/types/types.jl +++ b/src/types/types.jl @@ -47,6 +47,9 @@ abstract type CoordinateSystemTransformation end abstract type CovarianceEstimationScheme end +abstract type NoiseModel end + + type FundamentalMatrix <: ProjectiveEntity end @@ -77,8 +80,6 @@ end type HessianApproximation <: CovarianceEstimationScheme end - - type Taubin <: EstimationAlgorithm end @@ -94,3 +95,8 @@ end type CovarianceMatrices end + + +type GaussianNoise <: NoiseModel + +end diff --git a/test/perturb_tests.jl b/test/perturb_tests.jl new file mode 100644 index 0000000..9a58491 --- /dev/null +++ b/test/perturb_tests.jl @@ -0,0 +1,23 @@ +using MultipleViewGeometry, Base.Test +using MultipleViewGeometry.ModuleCostFunction +using MultipleViewGeometry.ModuleTypes +using MultipleViewGeometry.ModuleConstraints +using MultipleViewGeometry.ModuleConstruct +using MultipleViewGeometry.ModuleNoise +using BenchmarkTools, Compat +using StaticArrays + +# Fix random seed. +srand(1234) + +𝒳 = [Point3DH(x,y,z,1.0) + for x=-1:1:10 for y=-1:1:10 for z=-1:1:10] +π’Ÿ = perturb(GaussianNoise(), 1.0, tuple(𝒳) ) +𝒳ʹ = π’Ÿ[1] + +N = length(𝒳ʹ) +for n = 1:N + @test !isapprox(sum(abs.(𝒳[1][1:3]-𝒳ʹ[1][1:3])/4), 0.0; atol = 1e-12) + # No noise should have been added to the last coordinate. + @test isapprox(sum(abs.(𝒳[1][4]-𝒳ʹ[1][4])), 0.0; atol = 1e-12) +end diff --git a/test/runtests.jl b/test/runtests.jl index e952379..4c5705b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,3 +11,6 @@ using Base.Test @testset "Rotations Tests" begin include("rotations_tests.jl") end @testset "Cost Function Tests" begin include("cost_functions_tests.jl") end @testset "Fundamental Matrix Covariance Test" begin include("fundamental_matrix_covariance_test.jl") end +@testset "Satisfy Epipolar Constraints Test" begin include("satisfy_epipolar_constraints_tests.jl") end +@testset "Triangulate Test" begin include("triangulate_tests.jl") end +@testset "Noise Test" begin include("perturb_tests.jl") end diff --git a/test/satisfy_epipolar_constraints_tests.jl b/test/satisfy_epipolar_constraints_tests.jl new file mode 100644 index 0000000..347588d --- /dev/null +++ b/test/satisfy_epipolar_constraints_tests.jl @@ -0,0 +1,86 @@ +using MultipleViewGeometry, Base.Test +using MultipleViewGeometry.ModuleCostFunction +using MultipleViewGeometry.ModuleTypes +using MultipleViewGeometry.ModuleConstraints +using BenchmarkTools, Compat +using StaticArrays + +# Fix random seed. +srand(1234) + +# Test for cost functions. + +# Test cost function on Fundamental matrix estimation. + +𝒳 = [Point3DH(x,y,z,1.0) + for x=-1:0.5:10 for y=-1:0.5:10 for z=2:-0.1:1] + +# Intrinsic and extrinsic parameters of camera one. +πŠβ‚ = @SMatrix eye(3) +𝐑₁ = @SMatrix eye(3) +𝐭₁ = @SVector [0.0, 0.0, -10] + +# Intrinsic and extrinsic parameters of camera two. +πŠβ‚‚ = @SMatrix eye(3) +𝐑₂ = @SMatrix eye(3) #SMatrix{3,3,Float64,9}(rotxyz(pi/10,pi/10,pi/10)) +𝐭₂ = @SVector [10.0, 10.0, -10.0] + +# Camera projection matrices. +𝐏₁ = construct(ProjectionMatrix(),πŠβ‚,𝐑₁,𝐭₁) +𝐏₂ = construct(ProjectionMatrix(),πŠβ‚‚,𝐑₂,𝐭₂) + +# Set of corresponding points. +β„³ = project(Pinhole(),𝐏₁,𝒳) +β„³ΚΉ = project(Pinhole(),𝐏₂,𝒳) + +𝐅 = construct(FundamentalMatrix(),πŠβ‚,𝐑₁,𝐭₁,πŠβ‚‚,𝐑₂,𝐭₂) + +# Verify that the algorithm returns the correct answer when the +# constraint is already satisfied. +π’ͺ ,π’ͺΚΉ = satisfy(FundamentalMatrix(), EpipolarConstraint(), 𝐅, (β„³, β„³ΚΉ)) + +# Verify that the original corresponding points satisfy the epipolar constraint. +N = length(β„³) +for n = 1:N + 𝐦 = β„³[n] + 𝐦ʹ = β„³ΚΉ[n] + @test isapprox(𝐦'*𝐅*𝐦ʹ, 0.0; atol = 1e-14) +end + +# Verify that the 'corrected' points satisfy the epipolar constraint. +N = length(β„³) +for n = 1:N + 𝐦 = π’ͺ[n] + 𝐦ʹ = π’ͺΚΉ[n] + @test isapprox(𝐦'*𝐅*𝐦ʹ, 0.0; atol = 1e-14) +end + +# Perturb the original corresponding points slightly so that they no-longer +# satisfy the epipolar constraint. +N = length(β„³) +Οƒ = 1e-7 +for n = 1:N + β„³[n] = β„³[n] + SVector{3}(Οƒ * vcat(rand(2,1),0)) + β„³ΚΉ[n] = β„³ΚΉ[n] + SVector{3}(Οƒ * vcat(rand(2,1),0)) + 𝐦 = β„³[n] + 𝐦ʹ = β„³ΚΉ[n] + @test abs(𝐦'*𝐅*𝐦ʹ) > 1e-12 +end + + +# Verify that the algorithm returns the correct answer when applied +# to sets of correspondences that do not satisfy the epipolar constraint. +π’ͺ ,π’ͺΚΉ = satisfy(FundamentalMatrix(), EpipolarConstraint(), 𝐅, (β„³, β„³ΚΉ)) + +# Verify that the 'corrected' points satisfy the epipolar constraint. +N = length(β„³) +for n = 1:N + 𝐦 = π’ͺ[n] + 𝐦ʹ = π’ͺΚΉ[n] + @test isapprox(𝐦'*𝐅*𝐦ʹ, 0.0; atol = 1e-14) +end + +# Οƒ = 1e-7 +# SVector{3}(Οƒ * vcat(rand(2,1),0)) +# +# β„³[1] - π’ͺ[1] diff --git a/test/triangulate_tests.jl b/test/triangulate_tests.jl new file mode 100644 index 0000000..9a27898 --- /dev/null +++ b/test/triangulate_tests.jl @@ -0,0 +1,62 @@ +using MultipleViewGeometry, Base.Test +using MultipleViewGeometry.ModuleCostFunction +using MultipleViewGeometry.ModuleTypes +using MultipleViewGeometry.ModuleConstraints +using MultipleViewGeometry.ModuleConstruct +using BenchmarkTools, Compat +using StaticArrays + +# Fix random seed. +srand(1234) + +𝒳 = [Point3DH(x,y,z,1.0) + for x=-1:0.5:10 for y=-1:0.5:10 for z=2:-0.1:1] + +# Intrinsic and extrinsic parameters of camera one. +πŠβ‚ = @SMatrix eye(3) +𝐑₁ = @SMatrix eye(3) +𝐭₁ = @SVector [0.0, 0.0, -10] + +# Intrinsic and extrinsic parameters of camera two. +πŠβ‚‚ = @SMatrix eye(3) +𝐑₂ = @SMatrix eye(3) #SMatrix{3,3,Float64,9}(rotxyz(pi/10,pi/10,pi/10)) +𝐭₂ = @SVector [10.0, 10.0, -10.0] + +# Camera projection matrices. +𝐏₁ = construct(ProjectionMatrix(),πŠβ‚,𝐑₁,𝐭₁) +𝐏₂ = construct(ProjectionMatrix(),πŠβ‚‚,𝐑₂,𝐭₂) + +# Set of corresponding points. +β„³ = project(Pinhole(),𝐏₁,𝒳) +β„³ΚΉ = project(Pinhole(),𝐏₂,𝒳) + +𝒴 = triangulate(DirectLinearTransform(),𝐏₁,𝐏₂,(β„³,β„³ΚΉ)) + +# Triangulating with the same projection matrices that were used to construct +# (β„³,β„³ΚΉ) should yield the same 3D points as the original 𝒳. +N = length(𝒴) +for n = 1:N + @test isapprox(sum(abs.(𝒳[n]-𝒴[n])/4), 0.0; atol = 1e-12) +end + + +𝐅 = construct(FundamentalMatrix(),πŠβ‚,𝐑₁,𝐭₁,πŠβ‚‚,𝐑₂,𝐭₂) + +# To triangulate the corresponding points using the Fundamental matrix, we first +# have to factorise the Fundamental matrix into a pair of Camera matrices. Due +# to projective ambiguity, the camera matrices are not unique, and so the +# triangulated 3D points will most probably not match the original 3D points. +# However, when working with noiseless data, the projections of the triangulated +# points should satisfy the epipolar constraint. We can use this fact to +# validate that the triangulation is correctly implemented. +𝒴 = triangulate(DirectLinearTransform(),𝐅,(β„³,β„³ΚΉ)) + +𝐐₁, 𝐐₂ = construct(ProjectionMatrix(),𝐅) +π’ͺ = project(Pinhole(),𝐐₁,𝒴) +π’ͺΚΉ= project(Pinhole(),𝐐₂,𝒴) +N = length(π’ͺ) +for n = 1:N + 𝐦 = π’ͺ[n] + 𝐦ʹ = π’ͺΚΉ[n] + @test isapprox(𝐦'*𝐅*𝐦ʹ, 0.0; atol = 1e-14) +end From c470e8a2bff24750be79a1a0f0d8a40009b06e27 Mon Sep 17 00:00:00 2001 From: "Dr. Zygmunt L. Szpak" Date: Sun, 17 Jun 2018 23:08:30 +0930 Subject: [PATCH 2/2] WIP Triangulate --- data/teapot.mat | Bin 0 -> 41369 bytes demo/example_fundamental_matrix.jl | 4 +- demo/example_unrectified_triangulate.jl | 111 ++++++++++++++++++++ src/MultipleViewGeometry.jl | 24 ++++- src/carriers/carriers.jl | 4 - src/constraints/ModuleConstraints.jl | 7 ++ src/constraints/constraints.jl | 47 +++++++++ src/construct/ModuleConstruct.jl | 2 + src/construct/construct_essentialmatrix.jl | 11 ++ src/construct/construct_projectionmatrix.jl | 51 +++++++++ src/cost_function/cost_functions.jl | 18 ++-- src/draw/ModuleDraw.jl | 3 +- src/draw/draw.jl | 78 ++++++++++++++ src/estimate/estimate_twoview.jl | 2 +- src/noise/ModuleNoise.jl | 6 ++ src/noise/perturb.jl | 15 +++ src/triangulation/ModuleTriangulation.jl | 9 ++ src/triangulation/triangulate.jl | 56 ++++++++++ src/types/ModuleTypes.jl | 6 +- src/types/types.jl | 13 ++- test/construct_essentialmatrix_tests.jl | 28 +++++ test/perturb_tests.jl | 23 ++++ test/runtests.jl | 4 + test/satisfy_epipolar_constraints_tests.jl | 86 +++++++++++++++ test/triangulate_tests.jl | 62 +++++++++++ 25 files changed, 647 insertions(+), 23 deletions(-) create mode 100755 data/teapot.mat create mode 100644 demo/example_unrectified_triangulate.jl create mode 100644 src/constraints/ModuleConstraints.jl create mode 100644 src/constraints/constraints.jl create mode 100644 src/construct/construct_essentialmatrix.jl create mode 100644 src/noise/ModuleNoise.jl create mode 100644 src/noise/perturb.jl create mode 100644 src/triangulation/ModuleTriangulation.jl create mode 100644 src/triangulation/triangulate.jl create mode 100644 test/construct_essentialmatrix_tests.jl create mode 100644 test/perturb_tests.jl create mode 100644 test/satisfy_epipolar_constraints_tests.jl create mode 100644 test/triangulate_tests.jl diff --git a/data/teapot.mat b/data/teapot.mat new file mode 100755 index 0000000000000000000000000000000000000000..a0d30de98038a4530ff048dfd289646ae5bd251d GIT binary patch literal 41369 zcma%iX*kqh_;xii;{3rHA`wajlUFM6+QvTkdt?XroEOpK1P@$Ztv zYmxa+-A0onlv?y6#&c^&KK9JCe4p<9@7$!->VNOHP!r#JDp?r>OE7>R#qUE1!v->KS`;#Pv8I3-8r>_SbYmtm7R0>%-%-JKDc$sYtU3$gDeZ${03@I ztN66r=2Kt-?(PiI>@`lPNiHQ?i`IT+n^uDOpEv9~w)yBdZWf*!>GN?M3j8h^7f7mA zwSICo)d5p6CbRUouDp8lzDwQUd!j^f#o?YB-bbdhcFf!9c@^*=MFA;fn@Bt4zx-(* zt9QyCzA!A+C4T)lXYfd!zqGZGTG1DE*vMm)^OJE!Hxby-O6~XS0!Zhh+v)xLCd zvL;*L0F=!+`+h>^ES>PssAEn2)n`@fdVzF1jI3g~$9Y+&TS;A=scoTX%A~S5K2>?Y zc)gofo(Jy+zIaZHoAg+cj;lF}$$cZv6i3`=??%A?I$U0ZK+h`UM~gpj6|Zl6z!jb(99iHxQLw9j0{xPH*<*yo>H;aIiPnA zWl38b^_*@34+blrMz~@->}M;62a?|dsalW1f-c9Q)X|ZUY;vHR!HKbwiYU<$cZ3Ac zFbXfD_AB(?@6dWDviDtewyPE*!niu}A@gU&ok{4ejXW==R+Mm|rPwHT`xNJG1oh3R zFp22A9Q1Ym7~MhfQw^dxLF&LW@JutS)K4RPuK44c!S%X%#Km9I1=kX=CCPYuRPy5sEQ z-l$^l-@javX2p>n@RLY;bF0&c09YDUmdX)Ka?3_gb6C-N4_jJ~N?RAKtE1NNUyg4< zYb|d6(gE5aBlg=DQbuAjhs%t92cFFW$?0pFBv!_EQqSQCdI$Uy8>)zbrCUPT5ayN+ zN0+&=^Quh(K`MUtGKsKxFyDI9ce z+VVTI##P{OnlLRT1rd2$d`d>B{#<}+)w@u{BYZVk^1w_aOKp7_1~Siu_*ii6CCrs8 zFBbFU@v=vwE{T_x>3`SVoju%)X0KW5C$`j0gq9(Ep6Dxj!jri6fXrw0Cn@#Hi8TBx z{K=1yg=Mk$g3ETmUec%bk94A|xz{tbm(k1(*VOYJ4R*&_@wOH!@Dn4M@%xGi^Lw1E zyi*TmU`1+C%(oUArxB@hrw?0aUhklUq>Nf{n4sm8Oh&hpw6*`{L!sXY?cg@WmCBQQ zvZF1h;e%si5%86e!7|e~DbC`wb#*?md6*TuBtm}eMp=r?Y`1Wi|OtH-2aZ^d1N z!^((y#E#a-&OiC}lD)=hW0`FT=7?$AB`wB(dwHVdZv&ID`(60GyHOX`w401x+T)D+ z{4P6i2Z_^AHaCS>-I{^{wjgoVUl-QnTeyVhrm~CGbOBG+rI=3C&b_9q65ydjzlh+% zzA2hJE1&5c-crN#(F~(IPCA_cIdP+#fl3{Y2_D&fX~v_&3?ASM%#27ezq?x`1!-qV zW5R{qd@lmM^`Aq$5N;e0bH(L8u{|85lQ}+m>h_brXQb(ykyQirjApW^G&6y*O?f4^ ztU!-kk~@khOrkG9rQs(*-@lKrJ8;%-`^{+}%Z6K->9|GRQsf8i8K8KhblQP0i;q+G5X#jBqnub--xmcY_y&1^_3{1}4ZC&bFZ%{|_(kRW{X?@H(2=jPgpue-|FB9IWaG!#D2J$w@)E&o?Esct?D`B_*!@lCkn`=s}0{k z4o39r^RRVj)8WgnGG8F*v-ou_iRL`owZ~PziYhjm-3}bL=W$!CZKv1lP}?T8fVT>r zS9F4NE`pB$?~oeU9_>vhxc~1f?{%h4zr_celqNG|34@`X2VXp|Eyurmun1J;p8_C+ zC((+FOc&f^4pGy&#t#&Uu})f%VR93W!1u<3IlyoX;*Vp-&nACF#KSVWn1A^_Z|t5G zHHLgaI@5V}9s~;r=u^0Y!|GmQkw*+u1iNh_cAh0HO%XNQxLwB$p}!HWG~?BV1`A!` z3;Pwc7?B-pJoFo6IId0Y*0}5$eBKK7Kfts_xgZ!KPYjz|#<^oK9(kaw&Pgd<#f;UO z5bQSOBlC|q(jP3i$;GU|DbJLwcaP4DQ$_GrhL0t@d$iV|JHe4lX-_gvgVoM*jD&Zj zt}7QU8pRa%@9%wIZE(|kPqNNNodPuzJ`eCk_rOL7#tKd7;%2!n@$pQfeY|43CgNEe ztaZi34tV%gaS^zDIRq0ubF(F33j&wxS1zqybArp@-0toZB7TWuYQYEG9*0K1gEK+CDy7rmP?m0IF6 z06875#-Ynf;2A)POefaA`^$UZ06a3sWee*+>3fEq>0F}1qv8YSPZ>WQr(k27i zFB?s8DfrgsV#1KP~(4oDv!hv7p&cuXOfEByJ$DadzO&o{c`@615irvX7L1oN=9X-6tIS3@ib5w8nj~E@W5T4f;1;mHFIgp;ZS-Ox!N?9R2KJ-e{aIgKyKz zsS`d^WVRt>JJinlwjmTz2ff4kcu0s~WqZFCB35SxsdovlY}4N&YL+H1v_AwJvx^;h zH#D9Q66rikdAvjNM+u1Wh;;p?N%1jn`!2N_#1gi>Yq)OZXOyb-bYo;E(oww3tDOCM z#G;ih2y$BrK2y%VbN@KTJQk-N{4+vS`miV6v<$^F-@H^kfjHg2&e6D&&vS;=M%$Rq zXkZHmiFn+bhWDbnfK~Oq@8w_)R_+W)WCztvVZtqX609~g;>bC+BCZ2hHPwTP}R zjxoE+F7$%ULEb~+#|!M3N&6Z!!WWh}KYubKAH)`9N;we0^Om_bEpj^0;c$Fh4KNQ# zfM3@h*32FlE|HwPrPi_cu@Vc=!gg*Zj4BTaBdj>)=5H>K7!RPFUHI?zgUo(<*50HN z#KI8Ww~jV|E$)06$LX+?rMc&WNrHLvY#7me_rJrzx#zo%t5C1GUf`o{7tlQnBq_e` zT*{pHI^+rfsJaU&wAthacJi(E55#f#sA<0~c^k+XN_Yt*bbr~wwIK_XqGE!2tDAZ( zS^Moz%D)2k+buj zE0BI#E2DcycMPbE0)K8(ASw>BDj$@0rO|esMVf$|NHRopYv%*ijpsRwmE|KCMO*Y? zZ#5O3aH&&z{r7&dg4R1E?&|V!v|H`V!Q6Ar6dR_q(d`4dw!`j*H>>LDip$mA4bYh5 zZ45<}@++gMM`|pMPa#HszNt4!ZulFGpkVW+Y(;nU%k_JI3Yix3MqO%ks1}d-7bw{M z2zW2!6!YwZEGqS(PlGJU?e6)p`D)uS@2mC_8`7&6k7|fBC)r1j9AUGAEd1X;8NUqq zrGG(H-*i~~t&}Q&gI;RrbM!ppqiE6KCcH^HHX65EDTrl2O zezxus?5k5y_#vQeq9GI4zg_RrbJ!2^LpY-A4{cnwRT*P zn*RLa*S!Bc9nbC3=q0z7mz-VfKdk6ft7RQCo#{4c(VR^k%Rg04C3`pSw1|&s0M*h0 z`$(qQ#IQOj%&WHcUV6Pst)My%K))I6gS7D5qRu*;p9I9-e!XYD;WV%N{6(VHtmb>~ z7*fCfRRhbM@=BF2M7HIu(>L}_i>@TK(9*-raAsd}#hXQjC(0oLpr5>c8vfxA24_m{D7)zN_Qs5^>nG^b=`&PHua>fAB{?r)_XF)9szr=XZ@o(3B^(&j-<^91uoJZu2~r>G&Hm9qcVkX8(1(c;LQz ztQO&8A>IExLqq?0F37cuj|kh*7!jdE&@b*EK&fw;_1N&Ba_)Rq+V;T*8^KGC7Tz|@ z_7k401$*@Yphu60PH?XMTxH3m)u?)=#OZ;n zRujPAAO)c?p9o^ynID*g#mNMGWo>d4{$->uM4MU$o z<=L|6UxNR#=n@gL%PIZM;7`mnq~*fgGi+eAu70aV*6*NIc9Id{@YSV*6bs?g;297v zeW3KufVSTYC7}cR>fhVz`hC(wuB92@HYa&B&n_bPvGmw}%gwH4ksL=ps@uLE5#=fD zTRiGQ%GAo(<_eqXs2w(@;A)lKA~I_%dNrH~VOV8q8|%e(Mb=nDX8og?VT{i*ywbpQ z9`&p9u#DTf)@aH60a0R&U}1QAbu&OgoUUjYPjUNPa=v19l<=^NRd@(|U6`4ej%Kfx zb2>%r{MiqV>+i>gfC@Qu*Qwj05||0mF=wq^P)Sh6hv2vtX$#3TmC$ZPj+f?IqEHq> zx$4d=fRs6K9ofp@q%Cn}lD&t6`S#R7W4r-oW>Tl7hR;PEnMpIwV_jHNR}9&K=E*EZ zO_f~Q3LmAuDb=!`T&oqvo6*qJmY6nvrW*$(T<}6jjodW)&Yt*Q*wC8bVifO%0v@VR zUTlMa{rL)g%=N@Q{spP6ia3q~$xgM3-KZM#egEeBE}xgs7M1?Rq*ASEvAu4PTzh7z zo6t%2v$;M#N|w%5ywKcXp8cG-2vryMRC~?YjgP~n@y?Iyi_z;-uBu6BeIE+)O6F;e zMRvjPx2{72z4-T_OI%u&FPHC*L2sYed{c&u$joK@Y*Q?w7QH{us8?SA0_(31Of!FJ zf6V!5dov3t7lyZToQ=#VphLe-;84_6@ij}?jdQHZi3<{d|JsONgOFdS=?Z-iZL{In zSmAFI-J@Ca>gpiW!*qfluc7m!Grm;{@?YK}z!(-JQP(5fE|0YYr!D22#B6#eE#~x1 zH*q5gPp`5+@b#4Apz6tniyVHNU{1!R(Gbx0^5fSq#S3Gnf>a-bvTASz&-E~^eMAIv zl71}jo(=YQp6rTGve;pHV^gDZG5;jnsHw5#*%nHN7xx8O9geVs54bc+Fw!>3Vu+%W zrTUjS@h}VIJ6<5Z?e=fr>!hy&W#^8|Wpc_QhsMwq#=5EU4-R>k&SQafaK8%&+5Y`U z6+Z4OLdyAAcG8ZL*tA{l>#PldjPs4l`mOyVu}5d48KRKp}!B z0oLOME)43l7cO|R&)$r@eClYR0J2SFMbuqzqm^?Ge-1c-cQ*^}XR* zEOuvd0_SkY{j=wo6f3v?+EL7kHYv;e&#Ei81Y12+%|b(Q1lKWZVuDj<3$i{U>_pcd z=J;B7QNM6cXB^i^+AlI0<{^d3!v*G`JdCVFtMfrHn}GA|;eBx4ZgAtYW&`35bRu*F zG4PUni9*=4mn^1^Z&ZkDEdZKXIwMgg!8s(Yqws@*!yixT)KwrO_FyhclO8{kxxyL?3 z?XgZa;q-N(%WCSZqS`0-z!tNXZ5{=zl~Ulxj@<%>Z62DMw`rM=b)|f9G?XH@O7jmh zgeE@I5L)K~Dri5)t}O+YnDH*A8H*l;I9lYyFa^PC-&N5dNDSlE!oA;Yd2=Di$2bRS zpVTHy?gF#D`Md-Wk7b1vBc9Fi(?#R?NT|<1d>^aKsImFZh)5wv0+s=fZt}BElQ6nA z{qQeqqw&~?WJBc7ocrC z#MXAT(so*2wIHJ0V4;9xX%pesP+Xhe&Xnm$PeG$CGzf{mt>akUlKCn2WIq=>jFLyL zvgsq}wLIA@ZgWAjD z8$IzamE~wI%-LrOG|#_}LX%(92mRFU2E8(Pd8?gGv-!1)^^YpcOIJAlL*~%*0A8x1 z2mcSFB}JR5-!Ntut(^#7$MUxM$r3D{VfZ8~%2;Dw%YLE_FwH^i;)z}BEW0kfE_25V zF+f0?SvI?ZBoqYsNHCa09N9HcJkPVMJp{<3r}m1|p(G=XkE}=*Qt~Heb%uvrw2!*L z@!e5<6q`ZNgG!Zz=O_Op*as|auE8v*{QKJ-CP7I!-y3U^hH;vrbiV+|$@Y(rECfyl z^`1JPgEKNx;X$ePXQy}3B~ zyKJ2*Mf+C41C1Bh<6TfscItDlFtz~dI#DR6*WuJ4ju%GBxxQK*E? zJ+99VH-}Adc^5zmPHbeejlbhwjEL_spDMl-r|Bq<$6_!M0IA^496g2DfPG2G|> zkKy(TlP)Uhb0A9!kKa&nekf?5ld6|1#iJ}EXl(piDtqCAsf1M6WhLRt%3l#C&y(KS z+vk&aPmB8e?D#c!O=XMB;n&$LAJV%|P$ ztSdeI^1QY2!+LY}+WF(d?@3+B=SfxiIsd)>q-FhxA>9&SUO~DWZ`(S&weGw-ub-iR=rAhRPdI0?z241QQu+piD9j1=8AzJ`Wj(6bD|%5aJ(*om z6_+r*WYL% z&K-3BqM&~7tqK30b9v0W+Bb9BWyI{wdtPQEk^3B;wR-yso2$@(Rb!E{AHO2$uJ4Bj2cnR=q41~!JS4QrxPw~eo2QKAFUtg?2qXHD+gt7ufJ=fgx*#w z3O!Hzz0ii~z0Aj324-F!+ds|G0s#cFED4N=N9J3~kV3Tu#D1~UH_-m2t+bPcZLD)> zXfpsW6T84RUG0})_(5X>9qOwVGN0pH%pyd!ILWWyQ9U(RnID`Y*q@7~ie3CxB#BoB z!-;6Wd3%qY@RO3w28}nIAu0U7I5BNG89|oSdDEi7v&aES7m`G|QBH%t<=PX@?eKC3 z4juD3=9^4aSZZ2Wi9eG9i9xKJ%5Znq1@PSoL5?XachNX~E>xdAEYS5}eea3r-?LBF z#x*)b-!ng6DHV7JazVW%b^N7ks#9MN*}mT+T&qkna30?g+WS`zbJ!}ZjkUb7h+8pT zX1$n|^eQjSg)A?9zoY^e-VP5O9kFW@w*Y&U>FO2(-jwe&CJ0I7#@T=@HhW**V6=Yo zquKukHl-=5gcqadJky1>nlekzDTK2H0`*VX9^BStKwa+$kb4!=ML&2Ip66o^rd%Y4 z*pN>}yzp+M7DXVuzC8Lxax;`ga0UI`$=PI1B_4;eWIbSV9&8Rh*8b5p_(Z6X9^f1` z`a{{zuZ>&&dlJR5#{AWs!16K3Q1uSCG~<$^Hh75-Ydn@eI)TzzJC99hr*{~AgvJ+w z+Yr6-{z@;#r+ETKIufPHj^0$Ip;Ys9+`XUyA=hoa1H!AHI*OxIeujLer~%oCTo~mg zQ_$$nZ36$#ZRUc8b(2;E3-cp8fw+MMZew;{rY4HwJBW=Efm@rd>b6Gic3w% zQu8vuxI7-@+)Z%FA5EM#STIYY%8K{Q>cMCv*iX zF@K4a*Cjf}(kEy|Kys&Dm$`hJ45HCf4ez_zliU$Nz*|rVe9=6N`(R3x>j%J5&^*>< z*O5%g`}!{PH0=y19sMQIP7kAwSr}_jm%BjGG4{0MWRf4Gqd`=Wny4oh_lTjawWoXy z{LduLckez)a6!Iz_%*s2?^P(!XOu_}sgIBKT!SV8TTntg2NZU&&H>J>@NxNv$u7wT&ilM$>SJR-OB_c`$A2Rp2ULE4qC+_vaPckpwc5Iu`J#wdbe*A3Q~W-W7-(GnEJ zeM;mOwluJXeNeIzTw*+3mQ3v2n+gouxyYs78<|(YR)BL-2U}$3y}H~|AsLAa?2Z9a zQ$e@AkHjFh7|)>#SmZRTLCDIrAgvrs6Li}g+02$=|0W(?Wgi1gh+^jbmt{)E???cH z&i_%AAfnOqFWaW;6kM2$jH~M>*BppbcGmNRYC)szmv$p(gN!JS$?D2~spcF% zL|?ah-OJeCxldGrh%B6B-4qU88oh^3W*x^0h;NY1cO(({Jf_#c>%*7!T)I`(kbmU% zn23Kkb4CxWa6OFKMMBN3IoP63?h?7xS+bC*Xo__Lw^&p%-!8X&39!Z`cjkhJNjF44 zp{x=7P~5Bw+#h2T> zsgN-%my^t|XpnQC&~)K7Os(<}a!;CDI4ad}Vb?oYbSoue<~TE7;Xx&^8zeG zlnv7gx~FZm+xSx%Pt*%RoGAYpMsIdG628LW3_Ll^f&It`Fg{pdUHiezAaV@UWB%~3 z{SGc$d^~>ztajYPq=Wm7x7@OlFHWC*FA>@x=67V-#0Xr$w)Yl3PXDI~7_)mQrk25H zJuEB;GPK)?Q2{(m@_S7o{1lBZ>fwEM7qJkD`EKTQ{@};Os-Kie2gg`dSP|-r&!iOi znnlS+ge%k7VOQP+S1yDcdbjBTcMs-8B?M~N9?sBJ#{M$ZrgKEpP%)^!gl6ynw(;S- zT266#o7$gECrV0r@bSI$8OVQbjnaF2_|h8G2IcbIwb0=Kq=_sou0au>l9Lc=GqsG! z2X8XQ51VlolH8%Dhth+rk?o`N(?Q~f2h~#aQ8mUnUu;-Z3&x`4;zRBpK4$-iNB2a^ z$@6yu?;xnVw_R-5TYSrhc?W79w>Hc$Zd&ez={VI+>`6D@g!ab$@BH`S&%!j7o5qt_ zdOpGT7r7t{|8eX+p9USE7WMb)Z36!O;?_(x@LFU)bdn<`PD&z9tPWc_1$E0e`^(8) z_J48xfU`QULi(86hHB%-sQBWJq=-x&2}*aW_?BrQp|9lcv2i^N2m2F!x}(6(vf&ao zI!*<{FjJe&bW{1o zHaPyMeKa!-))@lnt3J9m<9bIO_mul|U@ahuFTz_r2pYoFroAmOty2gSap7sMgp}O@ zF=a}E1@_`R(}x+J@)tcO54fUGf0++hW0!Ylvk_+&i>*M{bCFp_Iitp*$XrTMnt9H% zH5@0si{%g>cNG4#LFT@1fu2Yz$u@v@zB0;$C&rODIMs_N2tnbZl zmNe(&M1liycCLg$NcX5EaRsQeCTKlt>l2$7p3UgrOZ6qFpBou|$TsuS*Gq|XQ`-(f- z{lKVCJ;`lzzlT)6uyq7bN{NC{;rak%w|VxTFmL-)k;1L;i$zFJtox$H;3VfZq&{;4 zK`rnK`RI59)e8nI%*mO@^XHiR^%911ML(JA4?z>C!EDt%6^(|LxO4>K+Fqua@?~#z zjFY%XsEmML1fsw}m!*ij=G9I!to{L=gnT*WIGFMVG?y9bzfRBG z9<)5T#6&aAPs5*@`|2_~H5C^$G^iLN?4OU4&5p3}58@ub3@VQpv3azun{-A#LH@H~ zOvH=dVn0*=9 z@1Zoaw|4Mw71@jL|66;Q?uFgc!slH((^ISxA*+}vzNz7)XtP0{{60o+_d;XiARml6 z1>*V}CBS7E7Y+r7iFGSQH!d3H4SdW}mS$U?QHs7P^a5uQCZdTleSMlbl5&M+g&v&=H6a#OmBg~t`CUU8n56bpy z$N6~zO5Sf38AL1w`E)Ea`_|Ex3gwC|si#6p-|B(ocMj-dm=Bij!2@~zZ+KWP>lJqu z!Q8c#FDu^gGR^+IGMH~qbhGdn;T9;F=&U7h5N)O3Z|jQBorh!sO@`yvbcN>H zzg)ByaHJ~#)cTvviQOOOr&F1~nSFvF84H(va|7Y6t)RBFF8n&^Y<<;35xZcu1~)$X z9@MPWtuFAX8Uyu7pD|ox<#D@T{@M;FO;=F|f3h{e!o^oj>!jy9GWS(RK%@Uhg#X6d zwUi}>coLv=IkC_h=`em`{WelJSlu_BP=VK}I4MnnY|+t-JW2F!jV z->hT3Q$1MRSm7?IzJkU|?wxE13mna)NFr}2Fo(9i=W_PG4XKKf-tl~anY*R7Lsx-! z>@9(8Xu{p)M~zv>{r};t6@4R)<(t*#!&!Nc^4YT02F>uEjlh?<6-Jvx;homfX+fY8 zTQl4DYxZYDDHe$(lLQ+jMs_J`#w9U&;N)R(V3T8hxUtYkX0N3FSrtfhX*4dPlr9=E z`=vWJ(W{DZ^)Sytc>?ukXYTkx<}4!`9%-OeqDuI7qVdHxuN~+sAmAU>DTfO6+F9q@ z_Rj5NoG5MuH6cfiw>mMsCuR4Bhg(@=)5X6-JFx8uVQloL^~-=o_7LO6$B@A1rW&-n z@B-F1XT`q1Fs?;3|vW-ul7xQ+wJ~maP)L4eC_Htt0H$S7qJ8p%x~kjtkJ1FG_2<1~9vEi3(FR-pQZ%~NRr;y~m+)*3)@=uANCz51~ z#wWx;*9#P|p9%?;hB$`U?X2ek6!&!HvJQO}TkGGySk_9#6`eQ?nSH6i=kMI+fzV*W zFzwo1I8TIaRsI>ob7op%JNQq6#RyzTqHK;<+ho8CAICTyWse@>z8|A*R_)wh8=ZNY z5`Q(Ri93L_Dw!W-+*=B*WfU&|eBt6-F692X>X%`5M!zre#O#bIhkyHR5s*{fJ`r4& zr2>Sn))f$R*JF$ z5*G9Jx9CE$8wPw-8mjH3%XU&c^D3hYeR6T`8KJ){d&>a+4XhUHvkzO*PGXb?}%+Zx*SBPBD}(7KR!U-1m6OE zj<9rDeSQKMHg_ON(2N9SF!igE_6SbI4a|ZB!aMN4Bc5H%*(6$=TwNsSE{b}t0R1pw z{GUpLwq5|eQApr+EoyrHko0a~g;sP>#0{tg8tMh@c42rihf$5+*}=8|>UYDIHwJ${ zXk8EqSJb&mE}?7wc?@qT;+CeSvF{uWzQj$sXMJm9M%eJ+kQmbVhctS@?}Jg)#TJ5S z95wuRZt|Kbh7vml7!}4eJmoz{cgZ~dAo!Lc>=?$3FE;NL`xDm)*k$afpP&IKhDM=S zJEPE|ORR)T^Z|pSB+&y#+s@jfkq4Y>_I(@*XkKQ3s{w1F@5XUd*KXFGYF(DfFU>t3 z?^wi_c-yQApcy`5PzxE~%=l-7Q@YGV*B(mmf_(HVfA|%Kt0OaKJ`SNlmv-^EOQc~< z#;q7oqbFal>cvq`&u-f0^FJ56d?X812!pg&XN7qjN+vrw77)MLHwb~*qJ=_^S%%K{ zlr^`#`|L+gbd64ncyvf+yoIkTiBcXjvt#oWdBF>+P@%I_yvixV!LL;|DYm|XxUpK|*U^b!sI z$kDwpUL~HbAlwiOhWb6aO0sTb3RCLN0zzu^WIA3=p&3Vj3W%)TiQmqzq@ z(aW@V#l+*R1wPgZ!|N%K#8F<7xRlO(kQ2B?Mk7m9@Tj%6dB1;xB;LTGX{XyjL$g>- zc_53}idpM2w z+Hhl%Fin{`adHuu>UCN3^U1Y6LS}EBiGh4oGakviE)%-+n(;&Q3@_{5|6}kzt+I=t`=oV~Q#3!!`sJ-R9S3J+@P`)H zH^x`EC=h?v>weN^JZWTq0GdfCue{gw=*{riUdU_q$oWQ2_E2f%5kM2O=W+OIb*duDcgdT3UVZYCJsR~?%YA+P_7RYuJk?HbC8DhtT7)mf%OQntmxH_ zRAFVRg#=GzN8XM?Ppg0KGM~8UpQNF9nFdn$wZCiT7OkP2le-3e&CWM5?)4e^)q-4R zCBeG%F)A%-K4mN;*L}f0m^idWG+^j@x@;#0+wI#AfJfFnF#>VPNOxX{{dNR8ypA1w zjlg$Bj`?ZrojSM$yZsifLGHdyyox#rY@r?&h94*fxQXS$)FW7tq}prl=_=46X|@7! z-GFME1w;@k3zhX-;Plc3VQNs-tE~_Kfw=h z{5{N}t{T|0JKgdDPsBefoeo|X8?7EChiiYGgcPl-(Z(1Yuxbtklz{7fxV@miBX$cZ z0ejwS=A|rz&B*La7e}mpsRSP4&}z~aSG#t@y|H;JbITqN8Gg!Od&u+!>^s>$!+ge+ zlz35D7Y|0|%C)u;{rV0s4sHhZ%z2Ab%sf1GAgzj826IDJiPAjOe+d*)N3ir>_~U9e z=5^x=>wdI(_JtEW6yyZ~Xqa-p)+^A>Yq0vTs)-%?rLDKC>c?v;XVAZv~oNr95T8ud;;}$Tv3emV#J(J3>ZPYb;a=S1Jucn z;ni|KD67d2B_OP8w7XXj1)@tDcT4|l+NO9iJHN;e7|iB#LS8=A+lOt4F1zeIdyT+U z<)RmXH~!im={ua1jB-FV16gvtZ6><9XphDxeDv+tux1i@t$o_4EiaqMd44k?px$Y^9@~t0K7$aHBDH7J zCaHu=Mu5`b>W_%qyZ|I)#*ofllefbV@b5u{qkbl*5pV7dB1CIG&~%8Mi2vt>jIi*t z;y1VqRk^`?L2oxq|4Q2_sq7Yuva;;#{8!pM^w=|R<-}g;x*sjL{4QTgVEpY~%+-md z_OIv2FZu;lzRLY_ijS4PA(-R^9QOm%-2LdCFt7#;*PEt>4#LqNe0seM6OwoY%e=~S zoTL^$xaoacZRx%A+7hQkAMjJoQ0P)byFt7h`@yGM24J)0wA zjpS;n_?36%SFO9PYd;Ynfc)6^)%$C`b3%-VZ+BjK*$62`cF6#Ajb-K7#3!o2FxWHI z`GMM;OJ;a*g6ltK^Ka}uJxp9hcAuZBSM?d%S8aRw^^L+1|9Jvg7GuwoQ4o}JQ{r)5 zDo$cBeOxktw=T88=fs$iyXD)&hoV|2aAO*T;#VpXAIz(1_G0Ky!zHIng!ix0lK*)E z8XjgLO3eCR7ZWah(mUA4jvBuQ$YaO6g4JQdwIJ2*OZYE`s~lkvH@c?;E{`K+VTP$? z8L}g)bC7q=D(yTx?E{|&Yf8~RCMas=jt*K0b+*l~JtgJk&{cfA%`YnxKY|z-^DpIsTI8o$SyYHopD(SMSmLVL>LwIxym49nkhk;jpL4f0SOTG=T$q933R+RSD37lg`o2o#WH& zAB^W+Y=QkeCCWKqq<794NY+`^MHu)#Zt`>V3Bb6NREQ(mnPrM3Z%@qTKRIh_L?A$I z!!Wk`r%)!UFVsd2(Epb@axRb!) zirILu2IhN{<~_P=rsG2!9LXw|o{dJ_8WR3u4IQG0+~OqXji(}S+qF>)rSSE1@_+l* zll#H4BAb8e?jG!~Yqgnc4F>R1PS6W!6ce`aub(;yKYRa~{-AGm55@_a_+9Z0{$)kg z?1=^FlNo`B*J`XNCvTMV(zs}f7fYbPB$2eq)Ax&91>V@b# z$#l~d8INyjqI)FnAY9w0etrJMz><*%9kf6om&IPPnUF6s4`WDQBBkmrEAGcP1$o@u zc~EuT+3CS}Lp=PpU(Ztt>R-&#;2R%hu{6GeHOP`;^+2`0+xBnTa@B}H8t4DbPU{^B zQD;BPCm)2>SEN|o!^il3;NAfi^JuRI?j)$-k=mdCT8iY{jdBPI4-%{IBb=dq*B*hh zk#N$Ma$Fg#gt5EYu&4c6ViLmY`iJI`^Q_&qi{+rJL8p-QV2?5o58|z&&vI0ECuRRe zGshR?o~dT1!Dfn|Bwejs8?LHJI0}lUJahiU9VAvVH#nI-4uI0{k_q-F>9!|R{Db6b z9y9%dt=mqnYhTAIh9Ta9-n74(Rd*;0q!x13LdLqT(Cp~d*FT!JRrjA0Lwc;6{YT** zWc<=-Z&DzY4{CC+Zi{)>Hou-o^Q>AIBOczDJS_5_jH;0!_(lm!0CVzOw4rQRGd97^ zah)2K=fT*n#R+*nzW$SOUU}Mn_NTw*&U-ltvf??)6_%vD@ZY78p|zu5y>t4u?|Wr3 z=i9e$Dfh>w!&ZF~6_1$&_ZpM^4Oar~`aDPn7Cej~W;A?7@eNNKpW?l;jJrIg|z;3aSLrM zV~apz3E3?_%u80TP1G90!=}c+B0V78u4tvzhWEfbXE7{F*-zWB9AhH5QKsSJE#)4+ zMf@}$SV9qBOiKNFj2W|gPcS*@SR-z_%fkvNKQqa8ccCB`zn!5E2FsYx)9WC@`|2Pu z*v(BjPsq>>$AQWIKhdz~dvSh40&uX#2z0C7&Xe5@6CXe?Yl4s9g-26nP<3OV8dGkO zj5mqTw{R7w788HZy?s4op*?i-H-Im?`hxu=IK=0n{zuZKs663nm7Z{rSvmLNEn!$O zh1373t)#=FI?Z3#qEMFin9oiE1m@SyQNA0DkjX2`oQSPyhK=G#yH+fzHnW0O4V(I9Qh?OW-7lF2qPKuh4v2Kll>uW7}Jn%oO^&c*=IRa7U zDt0g@XIAnRe0Uj)I44m{sumpwq@ry+I6v@fj{kWVI~w* zrm^8>%5q-=JXYZd>$9cfQ+AvnSDgurYca(wIDV=#+kK9I8m@VU!D1F5sw;$INNh{S z=#EE_59Q!Jk&9F-s`>?=cNuLAjb&za8_?HJ(iEQ4ES1&VgUUKCQSO6$GcvKIJVXtR zpiXYj1B5;!EdddO3AHin-4A&R+5l^Wy8Prlr9RqlSP+EXj0#u0B}9#^o#ib1S)^}y zaPsGL!to!Sw*LQzrY~`1`howh6iFpXD#sF%TnV}DE22m%Dmhn(5DU3yOYYo+B1S@p zIdZeP?>qN3w>dY*?B3_`dw$RJFYL3=`}2Oi-ml{?vPvs<6I&xg0F$mfS!xL-`vCHc zi$poDWc~zFy6|<`cmDk3l|4>s_9ICsV6TG1_J+2jldsN?$DZG|8&0YhxdiR)zruyF zM99j6d;F)eStHJ*P~lcu;@gwgGp^VM@t|-ckduHDQ`imrj&yY~?-g(rtED~L7*o-% z6&fIDg=Ma_^!nD!*&M0Odrvu`)CK0}qS(^%+?Y?{w?gLHnOC@zpIZQ3Y3;AA31Y{= zU&w`us{OkJ?>?0n-VNkxQrjkE<<>yUBGd8}TfjBXaW^=UbfOJgUYrfRp`UT>7jQp-E=%~iC%u+D!&9tW7owRZs^EK-{B->?vd%u!H#!(siC@`jdYv+99c;qo{(b~T{QLtCdHyT`&y z6mc8n2EJo2^$T3+3hUVY*KYJ1^?>o)c(;slSl zP(h=mXm!TZ>IV$Nj79%!MB2ieC$L}Iz8h_{dsZHgYK$W+DstqU=Q*W$3XpJi|fRw@VAnEml{VqhD)Vm9W zfz!MC*v;ecr{8fCNf$qJN73Dg>i2@Po?P6m*nRRwZ5fS+7Cxvujfe}dZ_)Z?ww47N z!bZpSd>hXY)?vIoMR#Nr`TU*3{ulc*z)(nI0XKDH`pGLO1Bw4B^MTFTYHHzu zl{tMsT`C$nKAF((Uc8{;cK!F0&zIs~`SNzZsw5g4RPe|FI$6De7JZ|nXapkmJy-{YiZ zqN2q|r?1s|VAD2R2>KM{pb9J+am`}%?rDajqZsn4_@?Of6c{ANqm?xubXaY=R0N*~ z@qo&oq}L!*%t%EqS!|Nk5T$oF(f>@ikL;k~g-LENqtH`)>YYU{jxjWA>(MXE9HMU& zx3sOiZXd%1P$=L{fb0KV)qKv6UjvzoJGut_Ot~z<=kiDU>1VqelAH+-e+OUCS&V!o zns(i;EaH7|brQbqm0^X5q-$Z@NUCvbagNe(az}AdQKjX+I(3Gi*|k4IT}pH4r#Gxw z_1A$OCHy@Uvh{qF5H%5V#3{ih?@DyQu|7L~7|Ju4N?&6w?4{aVj7Or&f_dUMML79! znhhrU>JCrk2hnw?P9Lno!Z$Z&rO2)F12uR00(^Y?+6hL*`AC-q3q}Qz!<-{n%Q3V&olxYf+ww?`EU*b*EziyqR zI60JIE0ys&TPT5>+T6wvKP`dgK=iBPlWZ8z4#Cf;2AS#La;v0_Rg1kL8<&M{bg{^4 zEym{Y#pNpS_}FFFYGJug2y=|^Y1U8j1tvar)2tEK%d{$q^|lnA1Rm!-i<8^Lq}3GP z^t^SQzP9XM2OA6g2jR+q93pq-&cP)>rxnMO_E48EGRLcy%~v>(rk8o()ZM73TLd+> zL$p<)*kPfx3J2XG`sI0c%M!@?#8+uz#r89PB3K5dwg5&iaD3|UUhoo|z`plA8E{4C zF1*R%F#1$5OHbm7@}BN(ie)+Y>mUKtiN*CDDC1Je=L&E)3niB7LUisPz;%c z4$dl`W|SUN4Z)+o|N0M|xP8HNxj@bJq;)Aq&6TJ>0_wSTbc}94`FvU?T5hYxq~LuJ zI(bV-?_fzt_1YekGvbqQM+h%8h^8YI)EM~l$ny0&FBM46MS1)&tb41;5uP0jiPSU(OAkCrpUhM zQ2#gmk&iKco^{cAZ-X)U%JJ>}Ov`HS2Ys1nlQ}7{`33-Drt#xK2@4G@fF05!tN#RQ};b>bf(`2QzUKzgnS%VIxIFd+UYNu($_&b04cu!oIA z_Wjeqck04m+HA3Iz1;NCwbv8OiL{HcjmQ2iDI7u8dP|m!61aeyG*g8<5k$9b)dE}$ zCx50Xrro{AZ@R&k#%)9UTD0!7jJL@9gn_bWfFlPwf%Z(x z(^cz$u*%LXFS}EK{l*q&Dm-;z7qSeYA^58}u>IN;oW!$V2>F?p`?XV8)vbr_2YH9n zuTQ?qv5iv{M8F;w3$Uq5x4dH>YJd)Jm`X2EpM`_B#)Kytj5Simk)nvxVhnLD zJ@{++f3LwxdutuwUBAqgioi^daxG0<;YwV`gjM|0G=@P~ZJZ&@pL^Gc{%c*V8A!fN zvU(4V< zg+aBSkUb=KZ_E-wUVCZZKXEY|c)$Xll7$9Q1?OY|OAeMjgCekTx-CDURbsMCI&8m< zt|Rj^Y78XBZLxF?$-1JyVli{|5cd3X$mfAc|F8=b0o-Ce=_|WT7K-B%s>r-^-+%rX zccZedF%uPY3@SYtLI0i!S^{myE0ZApO|D6c*PX!Ax&NWZ@jpu+Fzrzy(5r40h_ywlA$I6zY8>P$rYO2yo<)SCh7$ms<`NqhgB)zc&*X~ zcYxN`mx`|<*ef8V4tKm22p>TW9J^hvwLqvV*d}{#EmN7=+HI4uI3HD zT81#UGIune=m;h#d1fv&y@d{`3ac zQ`W9OqPI!~bQd%5ACUc|Yr%JM(kB)BCjdzj%+vz_eaEAj{QlVT*ZzvE-TYnSt)660 zJ)weo3!Ynsln`h<=o^AS2XXyk6R4e!ajlTCOrcdT_lfaNSnqJ;D!RyP z5bMS8r>II5K(-7G0@&+Dy`tLTIzJ%zCiFaGrBXDDqv8$UEOL*hvk+%?LQStscP(Yj z^jzRfQG@y?Ou(Idtev-9r^}4hn<3`~=-ygDdOW}SIa*isH|>MzO~EY#a`=e=$F?`d zpsrGI^W5@p_YtzvhAsKhvgQjy*S=-YdmAvV`dplTpT7{Pa02VBEsux(TD9dkhC$n& z9U1`EfdlrTd_lR3)Q;1%+fj|8hCa9*+5*Wvxr&%0E(G(UxKxQCFZj);vD|?IoW% zmCCNL&T$)s)idoOcMqUSw8&6IiK#EkS$lBl5EKqLf>nr^v6>eoIO}>}Y)`D5;yARG zqjx2044KG!noqvBfruDg+1F`ed5)jaLS0z}K?RV*`=^&pM2Q+KiGDyaOzpSnNU&6i zVJToEk(-Tu4uY}Q>oI3!*nM?2HLnKbj<^b?V`U~6f;{U4^%;B9_Aw>y!Aa_3_8CUopkW8zSxRy9GKQUmuP z|Ej6Mg2N9iQP&lf_V-pZ;4b7DuL_nQPl`2(2-qmLE&*38^PsU1&l|fYb$%e0Z4uOL zz3MNvcgiL`ubGgP7WR%wOTpckHytOpXAnB!25Np0UY^gdkmX7- zu({-3?5}EfKkzc;>l%4*w;pTxY%paDVBJT6<9Io+gC65N%(A6^$sapLBO>pkmkke6 z%7qYl5}M1Now<$jCxw56c)XQ*JEQPr&GzHzC9_peheoEdFB zbkN}L)rVU;vy~@5C}56*#nQ*XiNiIPP~mpV7qal!Ik^Y*-1M=6FPYw@NX68@y7j}U zn7ONd@Ttgta?W-=EWIRd4e7@6U%W*SC|EPxvb(vC%i$pT1SnK;pml<&pN%U)j4~2S z_ESoF_I~-lXN)<@Q^RHcQeehf3rsCQ{l{a#(oSug+@o zH95bu>6ZcX&-(vllEbo@uTciL&=Ol2a<6U~fJO#EZOyKt6dF>6AQX#f3~0F@#;zKz z8`8?%rQOs}Li0UVg>*w`4W42ZrPqG76bQ&~J2^P+N8e}Q+$3wDNw6Fuy)XthFKoUg zH&KyYk0SnQp#zIVpbLQhU!CJosE8%ROaLo4S*E+2T5A2^+m>NM&a`eqMcI?(-(pck zj}ylV#Vj898u1iI{LyKzpJeAZu_f*{UEi$lxM;tskSo6w-E`oE%>!ZA?h>d=-e9r@ z)}wzPK2Q#Ql%P1>3^{^iUp?sAmzKWs{+;o88FM$b8>Cs+?h{n})(aZ7w(|Rqvzoed zicNj_<(LH<%=O_^*5L=J(IJkQjx9yk$~ga&AHZR&yHksQpx7{WRhl9`2u=GCB`ar| zTq(53rJ)4#RZesb?0Le9n$z{)lAxLKZ#dGciUHE;D?r&9rHuZX^nvYA<(Fqn9}al; znmcXv)Nkvvfl^|4?3gY1x3wxrv3~4JTgfJvbH&H4X}Qu#?k<8W8>z&WH~x482Cg-^cI+Ly3Z>Z&}l4Mf;-2p)rd%um(!c91VJgpwZxr!mPYOvAqzfy{FZ5hVU<5)s7w5Ug#=za}oQ56Ae2E2TqKV*Q8M{xKUn7AhU>m}*5F zs|s9A*o2iIPw-qw>tM3SN{Oh~8}GC6Xm1gP)uivTyN!~)(VY^bW)b@npDSgwl{X!t z>t@|{G{&DbqHUBHAmE9&xBCNdzG8R>WgZ6U^3E#>p}N6Kd-Uy_jQt$)s*P3lPsj&M z9~=LDP>EvK1az8oU2woxA%PQU$kjUzrB9cBs}j)OF)i?XtwvPGs%IH!L7i*O{qJq* z?@W5!%BZ)$WCm55a6HA*hvvxM*bGRGX~d9C4^BSu#z;^)rMPJV1gBmxeE#wkptt2q zzt}zgwaRJ*oX=jgLMIcyuC(;gS&Mybi0Vq)1`7&NZj znTM4#rz}ktKiZ-44#~dm!6d4A+qn&-+WSZHc!R`Q-Y8#di347D#0Wz&4rZ}^J~NHE zSubzf2!t|;QhOO}BM5!*7R*<4hGz@h>A}UUfvxh32e*M-A;$H&R9Wo9(Aff|9*3a0 z+c}pPU|1LZrF5kl^~-?s_QNF|dW~zlhQJV?+x`4|(KB)u@yCGz)pMj1Pp(F%T!4?A z?2}(5VN%@WHZl^A+aZPH&%ue!@Np+~#J(m3P14-}{w2 z0(`9HhL30t%%gk)$RWXiQh&AbwyxU^O3$R>L1dfaNu0sQS2va!yRM3pNQdx2XqqKH z|L|Q>ic-42e(l1H_8g(rRDB3`YSp2=qKUS)&K1t6=lYO>c()oI(7t?uN-{e?$UnCI z*m>TvHlkJ8s^;3s+NDLUM)K=>Ltv)&DeC7{68-fcQqyQ#ZbgdNL#VuV0v^wL9{k7K9*5F-lv{-IEe*Ey~r5C(Z$b@+r_6LHGl@=xw42m!J zVmNJs6$cy#K{`O=niA+IxO=<5t1+iXK1(0VFkmkE`7g-bO!3oU3Ho0H*-Ud$_3duN zH1L6LkpGkOu~@rKd{c|8IXuaO^A|eL4)Q(48vA%?ED+9w3JN)BROAi(MxY@&)*buS z%~)qX(bLF{WZ+DIuko5K@tIa!J7MtH;ES~14?7+9Q0|tW!zEUzCV@yw;Hzeok<_xL zabD}IKUHRW2PW8x+auphH;Qy-lB6{1cp29!j2e&3or`nIKSU_yovtPNTI4B>dES@N zdss;!pOR;DC-N9io583$PFn9+rMBR3hh^TL%+%KcLs)37Ry4&S$1RZPuosrzh%5j# z;OgNE=QMTZADthQ*%L+b-^}Sbg+7>aS?ijVJ7iKWS#E zYz-@?1pNxae2wV`nrE3SJ#2gxLfz6&JCT)wxc{MB&L+?>X@x z3Z%%&U2#>@x{m{S6Zl>9JQL(M|C-`!z5U?A0?hyHo_mc)`G13byK>Hx&eZ5UKS|9P zBEOecglxmr0{ueOlc@9--GGC)S_i#R+!nl?z#loEW11!{-FPPJN2U$%lcOm!ZIKrY zK6Y5tv)zIVK|l8rAVruD3avd1gDFZ-&ccQVeNxglTL42y!$ZdN8yF#Ozwm51VppO5 zKE1C=!EL8(^#+n|xhKNN@e0w{^VW|GbIQB>2>#DiE)4m%mEZba5yBscQHU@hbZ z*sFHas2_Ic3RP?KT8x@W%`nZR8CXhc)F+*LfS|o(CLp;H&$dz1b3WasggE+mzi_Da z-LXr)Qx>$uYJOH)v}+XYnK?}*q7%%_{VkPXxB%@Vl~$tU?yM#?VlUeOZfklxerelr z5bS#px0Vr2=B+?!(8&{{{mI<3qGp31xw?TG_DQjS|s#RBhd@D+dx`x`Ovcx!rw47^1*Z=x>zu z-yD)}&^rc-?`7I4EklW47vMTq7zr1V!t@^Td=E)vY#gXdhTHQlemJB^7feDOlJBkh z5IjtcF03SlKsWNUZ6arh_5kq+8b zS#s+%wHM6vZVDGIIZzDKRkPWIWlfJaBgKyllY4bq50IgnrC~zfHX$x6FAWh_kNudL z5|1aTt5E`3vDP{wSV{2s))UAbJnaHEKBR_~@NYlw>PN?M{&|b_7f)6!J_`7ztO^_l zocd!)e@m;8ZDg?4M2EGo0Y1`M)Vjnyv^pohpFbS%^*ig^@E3k+&qA1fya?jAy2$z) zR$+5JOcf~Ud?B7^rnNp8c=k_WFLrl}lz_Y=MI!dI_#*L-Br_GRL$c6(8oSbKWW+3G=1zqdB>m$Cl_)?Qi9T9zBUyIx03`iF|7sCEmG-}i$@ zI!k}D>bFi=t}Jb6xMLa=h)Oi%eD<%ngN|=N`pQJrZQyG_S{#cT`}IV_aCu%!jGYB? z7;D@hzWl&l&niAFbx{g7_x>!UJR*t*at-y1P38l1M#LIl(_&-{a?&X#E6+_XXc4nl z8QQF6TzFbSBpbzHPq!_#Soy;xQq%gy1Bnq=@;gSarsSCExf&Jk>sjg{BIse;&vDSb zGv8NJ4tIVc;RH#VQ}^dcF-*5rNIcqN7-eQtmMr`m;xs^#f2(6qaVkwdWq7vvvK4 z?GCt?a(~3Z)R-;+nCM+Hfun_j8#WBDeetifh|}$dn)3&6QqA;-KzFIg96yEmj++KA zj@;-9;-G_;8dcnxn#W9ty@k+j-!N+IyrhToK}}|F0;g@7=3@sX#KF8?ElXlkxPo4s+YWA-;689qq+J7W z1RW$tABs>d#5s}?(xIwY<4CPN#WCj~t7lMW3%$#)IwbB%XBi}IfpFko%p?~xfN`e!8$o}6MwW$L zy14ZeX7ZR0<8xI>0RcQrZ4B*9G)z4;6Eh|;5dnN>D#Q0IwKwdWKE^=ZM%fICyJE*dg@dDE&5u0KZ=cRx*c) z87iKpc$XQ^o?$2BH*&+#LnjZ(fCXAitcO7W+^!`Dzpzm2=26 zyF-gKxBu_k?sA^hBVhJoij3QoAc|Uz3FjVObst*0vqFAf7h4irpBVmcmB5j(t+30q z&o0MjPbU8@&y~7|TVlN`=ck;{f(Jxy;4*Jagjp&-fQPMnB(BDuv%ZcrDU8)&B<0%G znf)6b4U?t{hu%~Fyo1D@qc&~QuK$aLLz8jpF$^=p@-ZNfcJ{lIk;6I2d>tE%{2X8U z2q6939OUjHKVoH=YAndP3oYKmO@IO(Sh?pgm`CLr7Z7SRfHw*$qCA5#PJ< zk8H>+kUNNf&~f$?T!O&M9I<-aM39k)H~m^_HujOVAV`6Y4V`lPyh0gb$Q1J|=9ZCr z?HB{K01N#8{b`o`qnAO7I)N!Sau(G%oxZ;fEc=ERJN@~g!B5Pk+8b$HXCB_Vt|s#P zzY}`j>oe58EQf|0$~yYK_#kCcZn>`duJ~&zD8NZXW6eL1tj^n@49pD7D3E9T`!(16 zljFrAUQ9T+t*`Iz+H>1}wX(e|c`>6dD4P7b<9mBueV527+p2F>rvohVENz95+ndia zQ?|`sKT1hm8@M~Cgw1($1-9&OG0Z2a*y3Z_aQ|m?yixQMQAXXLZxQ>{T=8q!u8HwR zV-_!h7c~FwSZ(3%Zt57;*8Y5eDSaXnplIBi(kHQ2rz$z%Gd|5;ZHvMy1Y@^$ur`dO4X6Er`s`2i|VrbFmwv@6X3Widq7k)SLCM;dae3e3L)U(F#@DHPrFK zZ+pvsdxtZ==!@g)_d~DNb5bY9XVtC|x-sYDAhNTMcJz1A51VPXGGaVbEFWKAlbh9K z-+XBqWc53wvyE?L1CsSK_nxk;C>Wbf^!YeV< z=EHSeFppBiJT8Y8@FV+4?pmzC_j8zO&>Ilb=X4WygIzfHf%oJ9=Mh**I#>CNu8$1f z517dT^_Ti9jlxDMumS{WpFfQ$3NNU^8HuMX6Dq~nB$h{q@7#~|=R543=8(z+`zP;_ zwK#nb%HT(5D=oymA!x}MQB#v%i6rH?3M}s zuAGy=jO_ev=-?{bpoK(kB2zQcx`mcw^_pXEU~apM1#^8_0^9xm_-NmUzAVuWK{O|1 znkr~@vz!;OTWOoO4`GM;SIkh`d@j9(@<$o$#^k^H$%mqv3~9jw*BYXFP9Fhufr^+q zfqi3`64mr>Pp&UFy|J&$7RS!)Jhcd7X9*CHA;6Yj9XShJ_s2KW;khwDw<|1(jS0jI znjCB6mT_-9*{NB5k=5ZPBTANvmi^Dt#Fy%WN55tY8B)!pIN)e>I+x)1z}5{$4*3w$ zTyH_TM07OK9{@(snLr-GKFS_g9449#LiV>QPoVPO7uUwj>&xp2VfP^cvAp!9e6mg! zKk(JJ-LK+=Sah%ZZ0<0EuI03!&&DW1~w_h0MDVc*26ZiroI6tEZVC-)R%| zOV$6L-bA${qfF7p!0D~qA7G0a$K$XcWg7r&z)8m72`I5eRfjuqe^cMIx!C7us|xDJ zsy-0CExN8 zZmY1TE;UF<{vZMz$=uh}UxSC2k_XFQj&?`TgpZ ze>>#OC1}do&)8Dink5lr|NWy?%0%FkVRx>}3H_Y1?{)zcIkyR%_@Yp!2-9FcGAt{o zv1&OX9DNgRhmdeX@BI=au263`f%1V6_JH{}b0`P1;TLtv!V0W@jy#D)9M;4uU znh{N6T;4>jdLP%tYN}mlDc~wQ(JL+u7b6%hT0}*B_tM=N^x z8Qyc3<1=h0^y3!NQ!oN>dVOy%PY?K}D*YX=Tz!T%bQsK`13oVIl%bZlq|sfcldwTD z+@2ogscwGaogD2O@KK~|z;g`}2kzQT|Ac&Rs$he4caz@W7k{cJoH-huxr_d*x7_C+ z`2o91`E5;o(#p1W1s24ORZp6VeZV@c`E5^s^(N6spyUdDKTG)$d_>OWJ~OkQ{5@EH zLpBLZRt+R%PjfP4s-G-Pdmrx1*S$Kacbw(?gjJ1f*i^JZN$;903z^4p zETZg|bQmTvpa@DKlp#n8Vzqun`6mBF45{-+YS_Qyi6%H!o!XJRt&J#P#$9@_?C&ab zXkuY5%KiG#mB7p5_7|}(Q7y_c^X|vSEumC;?P{E1#40#O;0sZD><@33;#@@ey*Hyo?15I@*CMur@uP4gSmwCJ@- zl&?~?(w6R{BL`h9>k-ZF2O%M-)3%Q_nRo5*5AI!%e>^$gO_vpHfhTlAt%T3A+$v6D zmnE7SIqV&jUlQ|3`PT^13~kk4&HcU9EbFu=@9pfe2Ea!2V|?a2Q#;LTWqfYeaL+bnWua@9sM^HZuy9x?9h$= zTc|FunB#*;dq_Z6Una59ORsL}Y1MsTUqCvGZYb2eJ+)(ct}KsAsJ=ABU&cD+cbCll zYmujf+htmp603?aROoLhKA>tB3;t6Czr!0U9p*n~N?xupns;p1dWKA;v86I!&xh{3 zIGXw1F@v<5xhMTh-@E=~*+hu{)$*7QuKn^tKCs6noqmy3sNw!$$B%B~$#2p{-75>^ zHuzjsO>T2*>7wl_{RUAAYX<8CSG1SNyLXFR77peI|cd5ALuqi@6(KMs_qT6$>MMW0t zN(klGss-J)As^X6x6%m^6Bp7=f46faM)23!&9IMTg{;=aE<#^Mb%F#FLxfZ_ZhWLl zw;(q1Tswgg`8r3W3HWsV7*~+4u$jY(0bX2~_YO)SkE4Ih%4k2f5-rnMp9tPKRy0um zs*GhVHJ?|ziu(xwA1<_w`>V%;3Il8had)Rd$P!Lkybk;+D3uvB`$c!=4Ufb;XqBFX zjOK2t?|IVc(9ZFeV&sOjHVaRjc?;?ThGXv|D$4Jg67zPo4R^#grg@<%0-Jri&(`@= zzU!bf24nAveOh#8w$xi*=;))C3)Tv-TICJ|&S@3)bi^~gA@>os^NI3kn3o7EhZO4) zMoru89Qga|nmKe-GhA_XvQPH0;|8OFxqyo%y|I+uMYd>_0?-MGBY+0*?S15cW^p({ zrJXtU(n48zsj3mJFvqEL4DP4JvwDoD&KeKlU!`t_GiV0T1ynA;K_fZixi+MZ;y4^> zK&u2+l6F~6{nr)73g(--x&y-76&Yu{6OpPuJj9`(I*YN~{(TcyAgn~fa`*kjz zOk^VP8&dGWm9(ch#38?e(x|+;n0>DdOPnsBl zAbKXG{)unT2$*x?3^s-BODolg)q#l$#C-rf!$GS6+n%Kr14Lmxfp+e0-@1m9;%x=; zH>8I+^6yR#0ExKVQ?x=oWb>8D*c~x26%;!UcIVyF@$!5CU(AnR27Y>dRNk^kk|muC zc;}xAP<;|u-~K!n-+ji)w^tm~EZ%+_wMvVK^$2QkRQA)e(L%dkjwih0+%sYRlfV;F z38t!q(~LVrPWJjlu z$RFRE`5zUtfxt;$?d{)3EB&zZf4grZD~`ZFcjbH*9j8`%d-iL$%oE9i7tOC^CPnNY z*ZhrsqjAl?a>jgW0&yOUZA+kwCzCeoLCiy@U_S`A=9w`-_RSbr+lR{%2G&l0zW(I# zcUexESZGC>)|cI#Q$GB`)5S92VPDC_il{$?tS(x;MuSt0q?K%aaxvBrS{ykJNVuAa zK3?A_uIHrH*ZUpb-NLRzn}RbmF0Dzjz;(6(16zV7wX2#ag)y~soUKZniMBPQy}`0QH`ms!E7&GLhV zlWeY_6+=7|cv^Ohu`)d{sxR_;Fl5jWCv+aO|32l-<79 zJWWtWM1MxKSe>>~_i(zgoEAT@4MUAh-BM%FpDYrYtCsbsH z^;as~9ZYIJD5{f-Nn5uA3b_`@_xsY$jHz2*bumb>R#k?B9o1_Nx-L zx!u=j`SrC?8*;FF%Tg|%S);3*t(b2}ZdqKb;^3(_{1Monlb0qUYt=qpWI2R$C*Y(f z(8;`)TCH`#caCl54CIQ&yjvOPv|2y7vR8<0ZA6YPc0XEXoa{NVMc?K#?JZ9UXXS8B z$0;ulF7{6#Gx>J8OI^(Ue-L&fxDOnNNcj>$T7u^Ni_lv@7*Du_Q{zbrI5STQjXO2#edr11J<1*E6SUcFV|Dd2f5vy`Z^*Qdo7+vlI}a441xVB4aSizY0-rl zsC8{Wh(@xXNPLOqMo2aIrmp72^(Lv{n*OD{@#EixR9@%zd$u)~eoA$Uu`8^K8$M5Y zaG{R}6-5nnu?Y2cbW(10?AvL-QHFR5qEsC$^HjB9r1+Yl4zOj9b#RLtxAgL@74&}D zC{5jBT<=6f*Ov?e?+OQ%chgt+mJ4LDQ2IfOYe*mOJZMqZQBoDG5Va&A#eEh2%l`sg zZ0(gO%aLT|2Hq(2L@YDP(o`UGS2YblNoSm$pp*L7QT26ZBjos(M!Dqut@tS*)L(fK zR`c#8BP=PM4I@95;)L@)GZ#DBz+KzHErxlxln_|eT31Q5rH*knqf~5c)ee!Nxa9>0 zNBEUtA&QA4Y?*@YP8)oH`*r-4w|}bKfo&(Y*+0JVl)EozK5lK#G)C%(@;}sNIK1rm zAA9GOu|n#&TkOCrPI)^EDkH&jC{s|Ws0EIc^80Zd6-pVVl^o%Mz+G0`xr?zNEvexw zeVu@oRcf0}dK+(sL|Pwv0x3vB6!;R;rcs>oLyTfIqkNV%X%3ZiJU#-p$Mo%bZph1R z5qy)EIj9*um)NPy8(Cqz{bGnn?s^}~sMIJm9m%OjOOO_OYvi9oebyy*^j&=31$Srv zo(rG}=a5$H_sKK!5|tr@mgd>t0n7$AR_}T8xqs{KBVT-%83zwajH zFA^pV!Vo_Uj~r5Y_l^k0A0D{CjB}(wjaEUH*xkZsgNYB+hOh@Ke)BSKD>aE}@|e4Y z{lC$dVtWMjlP;YVrBUrPV+)m@HS*{&x_^kxKsQB&$|T4V73BY(KPZldrrIY|hBTCKs>Q*I?Shz8+EUM^DoE!r!hVOUa)st@YT@ zYMU+Y&dsSN7KsuWV7(`Pi1wOn+}%!(Y^++$t2)b!)j1wFu>01c6I=IvqUyrsVE_2! zByBCs1GUGs`xgFNR) zTgISZ2X_#|&6y{xf_^WZ)9TC`%|L@4E(0mdWQI*0^$%L6ouc0vTk9^4%`&~LaF%(> zWvK42mj6aCi8PPP_N<*KWwlZ3qyx#SuAb3KIgx1B|7McQorXAkWWlZ#q0FvqbrEdA zTG}H-3IXyc@o|5=M6@_dCi=>8|ChxMzc#aZpb<511%ZmApR{agx{LYyrL7yGz=&gB zF6W(eG`DAYiJZSXYUz*5+)+c17nlKYd=*7|-hWzzC2q&J?^7?=DzxJA6%9Oa7=DUuBs_Yr#A3z*LyHLaD; zvcX+j$;0H?rzy9=`b1W4ADQbkFiz7J?h}t$>KONu6IsjCSlOi}zMnU>MLxgwH<*mu zIb<)sv#5~2DZe8(=k!VG72oPm~Un6+1p%|HoVj2>A;n{gmNA4ltkw zAES@>^MfV2rfX!t#CV4z6E+MUKT7Gaj4}DdHF~^(XS{7`Q~%p=$i~C@Cz99s$v63p zav5CK^$4_MtpD9r!kE6rW|YcLKBnGMUCybh0Yr3ThdUT*11iv=NOxA;xv+0u{%}*h zXMHg8JvMIXzaXaMeH+E)x?Hy(Hx-}KBPy3k{`JVO|Igf^gh1PO?6MNGFwe)v zYPEM@9dE4HhY~5+k9mHl@p)FCbY85aD-_>6NaoF|LJChojWu_)L0cIo!OCM{f5>{R zyhARhXpk=5fUky9|!*4*hGz{ z`w11qw$n_C{EsnaX#o#r0^$&;8ig>j6$;a~jTsaD+6mo`VucEttQG*VD|~}2U(!oT zgPkq!S=Dotr{6FCh-HG-h-zXY3G*CoZf6K6pPaKzj9Jl8){o2Q8-&h;rgLhV2pDaL zwWaF@NKMHuR+p3Q&2igPF_<5A0>{8Z8SEy2%FuCYbTp4@%;mgCxgTrR`46Il1D}K| zl$nX|g3ay*0p(S7M8=w_9bQiGu$|sDu8KhpKEz9prXoSp?&%%Ts4Lsx7<>@n{FHZ> zpu`noAcdtm2mal+fK%?43h45S%n}4!UF+VI5Jfmvoo8j4P1LJkv{v1WCYqLJJa9NK zL7lmacRzP*?kCVy+yV0t{#Y1zWL|J#oRJlrp>&p1;fk7<+Yrf*)#aZfn{dOzfNR?* zKEK`Cc_00G9b4d~V`7xgjZoUp;uQU)VMh7>MPlH`>NBJuVj0^mXuS|N2;;%2_b|sl z5k|9bW2)Z0H1!2#Z|Lj0BmP;`%gqk)QhFw!@uJoXo@<>__1yPHd~XcJCSBdJN)rg8 z&3(iEJ${ZT=QjV}{C{p((fru2EiW5UH1hEPJAJL2R110lH%-(aS9LLWmoB70;I_!Nx1rO13k@_2hod5QyG4<|8K45ioiLayeo z*26Q0x#(XH9!b!vy1A6shr52(%!@h*Sb>%%+XM?Gy^g~X&vM1BXid))I%q4uhEzIg z+>b2~9{qQzecDSRXdqnes1@uj(bz+p-++*|Y1Wymn*6cB^x)POPDWN3Uh9%c<3z`B z%&xbe-S5lj3YJxuytaIIn5_x*>jL$!&biB@_zt*0m}wgi6#*4+im4duS=Gj(wG1kL zN%1Hc-63zYjuQ zv;>=zkxLm8t6bi(!HfWiJr5N{Iro`i{7nv@^8U0H?VZc%_BmW{2R-JR5q?cbqnG3q zmT}Yq5zE6|A9bC)Y30Oo_n5o1+;v)f!}e=#RWW-~ecs1!Z4m-Zlbe&jC23%jv_a zbC&is)-~H#>|#iM!bHO*groO;P@Ti|<#MkeP$y4DDwX1wTh2rI*^hL@=a`8@Z_WG$eEwR(~I z>g~>*EiCdqr2Z5Trz~6&$`bGrvd2U_oLXj;C;DoC1EIhd2{40*fCD<~ia;s9g*JvG z_A(FkSJgt*mntx0TIEg{UVj@YBJgqTSCvsM%bwY$3vhCE2n<7G>+@=xhcUxTLc)Wv zPVj4xjKhU!P#rGybO=pEM8h^k)UnUMibP+yalkRXwPh$WNNeZ9&+4hWk2Mf%zFE9%QAkqHje*Vc)-R{_n)JSjgEZAc+^r< zg?=axON^g_wiaZsD^NXt32Iq#Vkd4R*`~qVJWlt3?xwPb!XR(NW)IZYUG$^-^p^i3 zGvr}};qOF{JO?7*cXZ9|WuyMoKP{rbi=KZ@iO__Xd4a6{@dZY@F2m$1rOECG_LUxk zXE*ZWI?rC+Us2nW2Oj{D2psJM`j&^a=GY09cqU7BDE8h?;lezU$g`|%-_Va=Xvh$eSdD;{Vom9p~@RI@B`;9;Ec>6J@s*M zKD*(uOLI~EDbIEszdJReu0V%Hcg6kAY~nj-CGMbxul6cyFi-pCv3@z%D>$ZI(x2CG zK*Y`C`(a$V5oMm)E~U(GkM5#HvnGG0?FzLT_HgUR6HOH9t^@ue=I;&99}G?0EVzGw z4Hh|gbh1Duh4y;z(j882a}&4Vm%lq2r3vsy{QTf1MX28#JGCl6vzKA{LDqsEpeP+Z zba2io_P3f25=gDw(#o??XNILmR$xYV%mXoa>rmQBG@{13NY@pi}we=NN zS2OtZU*Lwe&?~e$!W-|W+Ge`~-2$!?d#VvD)b)ua2!zE;8 z5gR6(zLcdG-^e*ew1wA*Qq98A1xq1Oh>@!+PbfNq+MXxvBSdNNt#i&M2wk~0HzVoV zBG7>Xno$17Gbgap@<~o9^PY{pGeo^MngI*FwwROWV!j^{h`2it-eP^-l>K}a(GIiN zx#<}bl5CL#sgUtwJ$G8#Z(0K_6?<<;u%CAe%5L&2$U`q>&OUU853TORa7S`a;UGQK zLt;R=#RNa)` z7xBrO!@E!-`~5)&4dH_1gtUoHO02F&)XNH;W!Nje*>6P|z%5Uv#>qH+ML2?tA}v3l zeYuJOQk6VlI;rs#;Cmr|TKnId!A~*?0xzcZ(z3V()6*e`eJaBy6QrWkN3)UVk$*YV zoQX`HHa90FjK;MUsQw((L^qK+z1O`4{kv2t?Eb0j9?wUDpnjW6nE_T~eYzK9!m7G) z)fiU3tPw5U#m@IPSzy;Ro!C$GfXMHx$(hkk&jfRgAH%6T8U2k7|AmIP&RX?6fVI zwYB2AY7Cw12~S%1g-mFgI$B7wdVy_Lc47kwU$R90J)A#;JVjHh(k!a}>U#Cipy*O9 zd{c2T_ZEJ&8)ImDd;k@U{hGlD=cSgfq3CnRN7cjogA$UxIQ!3AMRL0-Ku@ed@X?E( zRzW`&S@J^W_%Qj-e0ue7F}xAYNv4NmfOstA6W#bCuw?>3Z^Gt|0`?k{TR^y-G&Y(< z&pqJZ`hc52OdzeGy019OQxzf*!PTl8>@Niwz&dXqowg|%q!ATTC(SN#pz|JBZ>lJ= zgjmLWvn46Sj#iQFlF2}QCI*y(qr*xc3l9-uFG)>k(Bde@vDBhRv7hbwfWU+c`=xfR z0vrvvZQoK0#AOPK6$e!RdwFxfT`9lH;#y7&AP1t27mE)C2%r-5tFatbtj`~B0>f@Bj zr+b`rxH2T}1RJ29-mb@m!)|LIuI@d26@!}fGrB&7n8F9H!8W;6%&~@y1BBv>=*VfS zW{+LjE}$aKg@vT_$j5%m^|4uG|7aYx}9O9R^FsDt97`|;DP9KbWLMaBhUzibQmi!EPi1$aG z)Ou7Q|4DDGV&)lW1tq0Djh{W&+@Gi|!p`#(MfwyhO6X~;z1=<;h zW$lkUVDF~$Cb-&Mnx%VA_NUx8<-ch24T>Xfs zXhc^K#dV4~9#>VUAtM`OEra44O2I3~n8Va{Q9^(K7T_&2E}upTgCdKleL=G5e$H_` zw&Mu1I>%%!eRNdl?(RK*#6EsRM<(_gU)$4kB@2f1-{`q)Z}(6Hsn5e1pcSdv6Sps_ ziRfgFR`r|*T}roUC){DmpNS}H%N0;uDm1NUB_PcFJs%te|;*Z`I^c zXt;$VSvtELP|VpKreNV&VUGFl;|O_o^IZa2QIWH#ym*$|Ozc|OS6=#WuTR=NsCtqo z`^x2I%2H)aF98UHz+e0(`n4oy$qjS(|DRU!A119F^{-gGf4hEEv?%(i*X*dy$?Q^e=;tSCN^0^7b(o8_6PEI7$Rzw*^gBj3gq}`Q^x&$bV<4-1qCgn5wN#t*eydO8O}3H3-{W zBTXz?!hdfsqd{uD*k$U`Ze6I6*r2H>Not)g(N5*dC#+YR|5G9l<;(f{xSB+K z$ol>xVGH6oDo83r-#^6hMqW{*2UT1CK#UK+U-`mEHCJF@Z)x%G8O^?n`eF97LEcf> zy11RPIEX>ktY6DCm$Vr_O@^13;3J(5FH}>OX}JDLVLsgtA&L~2pU3~eB$pnNW7u3K zyDdZ4Wq|W5B7e(x4~tO~riVF6p3urqyLy39tCj&l)u;(CO%D&H1?p+Ko*jg!-W1#< zgraaJC{D`5Fy7TiO;<0hYRK5XqA4KTg~%Ex#l3P^?rO}uY^PB8zzN8YIv| zpjlYF8}*PH!3qAj6kcb$8Q2LyP~?n%y;*-_I0%}4yr=uq_awDeON3)dnk8*gHH>w) zS0%U&T53fxv>>+9``|mnjQ|X2lHggc3{u;51Y08L!!ATbuK;)!qt?HdNs?LF> zj%g%Lc5oQJGX=CiS1HQ6Kcljo$h}oses)r%K(8dnMCx$(7C--}!yf5`c<$@M07O;v zvu$PWc_FjLc0~uuI=R!+c*gVq>|e+WihRI2n^je(r}laaIVwi9r8Pz>k{g-=-^x2e z#!ToSbRG|FWlJ#lQ}SEu>j+0_Wvn-O>{<0FLOxZfEGxZ=;8p^WgKSB$4Io zC^3LMd!Qbk6U8xXRa|0fr|eXp)k>p3emiXYoYg7xd4XB%lEAt0*8L61+Q}f&WhBeB z@Yn25!Az{uDOy$tbW)r_-IisSH%uinE4LHwIN<`S-L-~|&Y5-5e^>fe&*t1CD}B|$ zrMGyN%&$a*TUggPDJM%EV~=Xj7byunxif*zVa7^DA-yEm*c;R*c<}-CYbQu$A57TS z+Fk$&?tdkA-BXcriV}?m-m@o&sp^ADw2Gaftv^~;Y1O8YV=>`Nupq8qKki!95L3Xc*1;P7pU}D3`}+P99}~K_)*3drS2?C z^%o?o8}HZ^++KS&al1N#d~$n5l-(l09XrVIYGVcq<){>n&UHNiAtAD8dd7gje5fl_ zjC}L53;tLsQH_8^{D{brM{Er2oJQL0pNoZUL(qRBR18fb#gagaU+Y6pED8kZF0XQb zY(ac7YH;m$X_TE<3(nQQ#OteJ3q5=Air=P@7CGDv2i#IaSa|h?m=?HAGHWmGPO-go zR)uJ$Q!4LMajHZZ8CoHEhjVZ6XqWlth(i-edrKGK{;M8~km>gy@$fhgC7NENG}}d! z6_O-hpf%P>=)ybWcna0?*1#gGVx!y8GWX#oR`*j`(&P?swGbUS~V(W@^F* zicC5=4|x6TZ9_%q&vFB}C`v;ksSYIH+$p=I%WR5?hzR4grCTzmA z=k$2Z@{zNDhLJ;a1O%9fZSI1bd;7L%*zq zK_4)jtz(gY8XO2NX}?1L>j1ks#=r>US~#>7_9(Ifrb@skQeFN)yX&VCloxmNe7Cp; z?5*6mFZtQ79g|o)RKnx3C~<13KBkhG*$RG5Eg$!_<=AiuLakZEE9456H62Yo`v!Z0 zbD!xBjyQkv;?%KGy%yzH@h9deXL=>3^LE`Mvy$m$gNU~`!3l`(z{TeF?lsUvwOk~CJwc1PDSBL_>;f7$%Oe&06z zl!=y1Pjih7R`$cFo%G@Kwm|_uV=Ts;&e>~~PlW$Eh}P+j6vv(LW}dNkd=7(G!xy3G z()|a`TUUBpDPY{W&2V|(*)aR&X_4ShEK9uEnmV#C)E0ikZ!#Qoch)XZ8TPPv=La3T zv?mUv|L13GKAT4Eo7fUi+B0sx|JuRcrU>o44f3=EDCLv7>H9T1TI3;%VysvU28~&? z?`w4H0Qcebdn>@~#0WR-V!I4d>ic$z^V4jyA>*_cc8E~CbNc;3SSavj(&9+%rZ8tR zKB`*R3*l;_yOv$R3Z5-F(3o{X^nFs6RcJ?fB}wihwPZ%V>oLUFjdK`Fk^#?KMQdf) z&`OBR1VXr$ekfhi4`cuuXUFaQea;H(%ywzJ+a*ei)_Lx@`%U`5_2>OwIcB`utEwhf zoDA)+mfmyNcYbE>oSh<(G^qNM?Ok)2`XEU8HuFZ6;AeJj9$(m>ub(qlPQOzq@M2Vgy5mMqU0>!0T15S&|`f0hINeT{xycUT&QAB37p1G<9lCV(ULGmHqPI z@2TYGWTQ2dM99Gildc(gupXW`Q*Bh_Q@;rL!M4sfooHqXG8Ng6?-Ht=m6)G*o@Gs2 zM*5H1PS^E^KT8{@4%{m3rP9OIx6ofc2g-X1F>SVv@jtZMK0K;-GVl=MAI#~S*4MtU zJV@JHzszOx9KYd3h>*-w!!Jo}%u&u67ZaAtpY%KwXfgili@Siob`TGmJqF;~{7|ri zX^~2Hr*_@51vhLdZe#sX8#m5`4^4F++C>ItjaF~XI^#6hBYM;yb(|9DdV1)f@y*w> zI-jGvQR-7_(#y4Qz=FcF?*iLaX52?mzi%R6wN?mhvvDQpuq5QDjA-r_DIRIZvJME@ z_N+Qpk&9!bs)N2jp~my4!ZZ|4`cENOO~kPe<7XrLPOq6ZPcA$#pxo~Wn{Ma9U58ED zCULrN7Cx^6YhUPmSjjHw0mi{5=8!~&z8!~cA%uv|4nBnQ?BTVdwS}3)y470XM4ZN?9~d!y(X!J3B$IK^>kHjpic1 z!YYdI|5DeRSk_e&GB6)!&n1b8q8M8JdhgGmY{<8rY?6#+L7qza+>i%=!s^-P{boHe zmjkn^Kp-X5B{zj2`tZ6%*Hk0}_NXE&+xJM5lb%CWe{;~) z#IdrKEEN^fE#_^=wjwuT#(oP9*x3}PTi;9Me2sEKY@I_A?8olgF{8xt1F1 zX1988b?vbn$VaoWpkD>ikh5TL;*cB>0Y9V=19`Lx5b|{ z;lm|w%d|bmC2I3??^ZQA+~*G(%l5oJ;G&cwQF+rP*u@b`S)X9R(%KNrwrT=}#L$6* z(&LvjJx62ooQZ2Qs`$6ypSJ1{Y-6x7 zq1MQ>&gqbOEUMU0|@*eCF=8#?acHySOq0YE-J~{ai?Z5^dg zB5b1>5D&{JyD%ew5oKo_KyO3&yoxMvT&Z=uM!z7!cCCE#k-jP;_7M+?Zpybcj(!bn8(%MNoHOA+iz7w;K-RX~Pq0qA=!OuCa68 zzm`~d{0X{Vfb!D$i5WFHebA&{K`&I-|K#8wA!D`3F^fFY-Irdu?RqK>6W1H}zC~V% z;?UA~@zYZ4gTZYvXfZj2iPyPiamOoeH4(COyX#SC>NBw|Yw2N$Jw}`b8^n+M$In(5 zV*&xP2iALTlZ#A%BsE@CEea#@M zAsJ1OuArG?63cN*pG0X?UR5r^)Wx9#B@6m3b#?VD>}kUl$mH!L7$_MxvqFpZSq%jB zhW+?}7oWOxlqfsXC^1>Z#@*Xe{`4CC+Ud6s{HRtSp(S3_S!Os!V+nO(gs1x2M@arm`q@zN+Q+n^H@g9<|L&Y4Xjib#c-yG>g ztVvS=5cNJlq*bx~4P(^%(gh_M1NNiV>2FV)@wc|+10qWo9?FsH8MaH-7mbA3D$R1p zFYePSPdQ7R{hDbV?&^)VYKGaDv;al2nRayHS$4Mg@7&=L4;WSM?0!`KyVMZN z=mHP6g*zbCdb!j78al=>apM6yV?0V*kY=K}pks6xgysq0;L+qGpmnr9 zHMbF0pNXM3RxNYhjsaW%?{+(-kp0Y)teCSAu9vr(2``;NbwO7bAyj&t4ajPIg#8(Q z5&yVfJ7t0lHQh&_6XKjMA$Zdu=Nxg1c_Eq+LCeDl&|A7AVFKAGKUEVE* zp`SI_ln;?T9Ry|%2>K@n-PVbM6h_01{@OPRb1GgM?`7@nNKx|xQ7>-J(8{J_QT~e; znGsVrt-AnBBGCV945=<50G{2^J%y6K6}6x| z|KAO^Qm6;Z4q2rBVl|%d#EcpuXv9(8DZh=qsu3#8-uR)3*hOn5GOylUJ;gZMKdW#4 zAi&`m`(<)j-f!BMGa2Aee}yBQk^QpIfB1J|0bk_RJlSVq)qc$VXsD))h1LjSE}yaL4!K2;xK0vcSL)A(VBJR1WTa;vE}qD z&t!=)Zn-_%EV-9IeCtx$F7QDd)(`k{2>FFFJ;GjFKZ%F>&Dr4jWkXfAvEC_PAXYH& zsC&~3`Epc}``b(GA#7s7Udo*|(_opr$iJxpkmKcOW*e-x!Qd5OMeZZ>!UC{#Uf{Kp z!O3~sfA#0tzk5807a}ef0e2pju>$@ntF#4|IkSgX#r;Oorz0nxAWDDgJ42SEtGfwZ z^r-KKw!mS;P430|C0pwIXkvRr@a;K#cR{^d^G1O=-SFA2?WpU#zZ@O7oPQ*b!>i%B zPxskTTqCaeeA9A(TYdeTg>u+L-AgiMpSoHwwiiRXO^t7B3s48{?hF?mOk&?56|uaC z@-IjBjW}@sKH-1JzA$?(n-%tM60#kjddweHC6uw0WjflnkO;WeXhH2eP6cjzvQ)I< zLals4k9d)dG0ex}NxxKF?RtKD@geL}?=sqF;+Y6M3KYzUXZx^FqUq*R7-XG|0r79$qi( zDOeqCzcwi5tTeIz_&JoVqbjGjBp?$Rca@_+4x12C1uMTS_{qMl0Xj~y|4n#ch6P;y zVVoHf+B23#M7i@t1FL8!iD-T)Lz$sfMAeEtHTCSZqkg@uOSc#8Q@MH8PlfdzA1*us zThKjAETIZLH%8rq9=^|YLX>im4ic;fW~m$TRI3*+8g2?engoM!%2{Axh7R4ZqB+g@ z8n$=Nr|QD47s6T>lVdKvi0Z=}E9$2MdMA`~2LW^9URZFaowLvB822H1EAb3HNT z^Tka^+Poy@?duzP*wI76rZ|Y^!j(sEvO-4%Rkji$Cm3EqJxjuL4;Dqkr(1QhdYyqA zqc5-c&Z7(}1eQ1-Va6jc6qjbXL9E7TB&Danp*h^2EE%he%6Yf+r@^XZU|}yK2B*{z zYW7r4apTawEpGf?=~l8;cs z@^-LIBtsQHG6q=hba=^LNKJ*o4zEK-Q{?avYxauFm10_hZ|ykhf45#ODC)W>!(r_e VvjYGvX|+D+ 0 + end + + β„³β‚‚ = map(𝒳₂) do 𝐗 + 𝐦 = 𝐏₂₂ * 𝐗 + 𝐦[3] > 0 + end + + ℳ₃ = map(𝒳₃) do 𝐗 + 𝐦 = 𝐏₂₃ * 𝐗 + 𝐦[3] > 0 + end + + β„³β‚„ = map(𝒳₄) do 𝐗 + 𝐦 = 𝐏₂₄ * 𝐗 + 𝐦[3] > 0 + end + + total, index = findmax((sum(ℳ₁), sum(β„³β‚‚), sum(ℳ₃), sum(β„³β‚„))) + + if index == 1 + return 𝐏₁, 𝐏₂₁ + elseif index == 2 + return 𝐏₁, 𝐏₂₂ + elseif index == 3 + return 𝐏₁, 𝐏₂₃ + else + return 𝐏₁, 𝐏₂₄ + end +end diff --git a/src/cost_function/cost_functions.jl b/src/cost_function/cost_functions.jl index 91ba160..278e796 100644 --- a/src/cost_function/cost_functions.jl +++ b/src/cost_function/cost_functions.jl @@ -1,8 +1,8 @@ function cost(c::CostFunction, entity::FundamentalMatrix, 𝛉::AbstractArray, π’ž::Tuple{AbstractArray, Vararg{AbstractArray}}, π’Ÿ::Tuple{AbstractArray, Vararg{AbstractArray}}) - β„³, β„³ΚΉ = collect(π’Ÿ) - Λ₁, Ξ›β‚‚ = collect(π’ž) + β„³, β„³ΚΉ = π’Ÿ + Λ₁, Ξ›β‚‚ = π’ž Jβ‚β‚˜β‚— = 0.0 N = length(π’Ÿ[1]) πš²β‚™ = @MMatrix zeros(4,4) @@ -130,7 +130,7 @@ function covariance_matrix(c::CostFunction, s::CanonicalApproximation, entity::F 𝚲 = _covariance_matrix(AML(),FundamentalMatrix(), 𝛉₁, (Λ₁,Λ₁ʹ), (π’ͺ , π’ͺΚΉ)) 𝛉₁ = 𝛉₁ / norm(𝛉₁) - + # Derivative of the determinant of 𝚯 = reshape(𝛉₁,(3,3)). φ₁ = 𝛉₁[5]*𝛉₁[9] - 𝛉₁[8]*𝛉₁[6] Ο†β‚‚ = -(𝛉₁[4]*𝛉₁[5] - 𝛉₁[7]*𝛉₁[6]) @@ -158,8 +158,8 @@ end function _covariance_matrix(c::CostFunction, entity::FundamentalMatrix, 𝛉::AbstractArray, π’ž::Tuple{AbstractArray, Vararg{AbstractArray}}, π’Ÿ::Tuple{AbstractArray, Vararg{AbstractArray}}) 𝛉 = 𝛉 / norm(𝛉) - β„³, β„³ΚΉ = collect(π’Ÿ) - Λ₁, Ξ›β‚‚ = collect(π’ž) + β„³, β„³ΚΉ = π’Ÿ + Λ₁, Ξ›β‚‚ = π’ž N = length(π’Ÿ[1]) πš²β‚™ = @MMatrix zeros(4,4) πžβ‚ = @SMatrix [1.0; 0.0; 0.0] @@ -233,8 +233,8 @@ function _X(c::CostFunction, entity::ProjectiveEntity, 𝛉::AbstractArray,π’ž: 𝐍 = fill(0.0,(l,l)) 𝐌 = fill(0.0,(l,l)) N = length(π’Ÿ[1]) - β„³, β„³ΚΉ = collect(π’Ÿ) - Λ₁, Ξ›β‚‚ = collect(π’ž) + β„³, β„³ΚΉ = π’Ÿ + Λ₁, Ξ›β‚‚ = π’ž πš²β‚™ = @MMatrix zeros(4,4) πžβ‚ = @SMatrix [1.0; 0.0; 0.0] πžβ‚‚ = @SMatrix [0.0; 1.0; 0.0] @@ -291,8 +291,8 @@ function T(c::CostFunction, entity::ProjectiveEntity, 𝛉::AbstractArray, π’ž: 𝐌 = fill(0.0,(l,l)) 𝐓 = fill(0.0,(l,l)) N = length(π’Ÿ[1]) - β„³, β„³ΚΉ = collect(π’Ÿ) - Λ₁, Ξ›β‚‚ = collect(π’ž) + β„³, β„³ΚΉ = π’Ÿ + Λ₁, Ξ›β‚‚ = π’ž πš²β‚™ = @MMatrix zeros(4,4) πžβ‚ = @SMatrix [1.0; 0.0; 0.0] πžβ‚‚ = @SMatrix [0.0; 1.0; 0.0] diff --git a/src/draw/ModuleDraw.jl b/src/draw/ModuleDraw.jl index a1c3046..ba3ac43 100644 --- a/src/draw/ModuleDraw.jl +++ b/src/draw/ModuleDraw.jl @@ -4,6 +4,7 @@ using StaticArrays using Plots, Plotly using Juno -export draw!, EpipolarLineGraphic +export draw!, EpipolarLineGraphic, LineSegment3D, PlaneSegment3D, Camera3D +export WorldCoordinateSystem3D include("draw.jl") end diff --git a/src/draw/draw.jl b/src/draw/draw.jl index eadf1d7..094e88a 100644 --- a/src/draw/draw.jl +++ b/src/draw/draw.jl @@ -3,6 +3,18 @@ abstract type GraphicEntity end type EpipolarLineGraphic <: GraphicEntity end +type LineSegment3D <: GraphicEntity +end + +type PlaneSegment3D <: GraphicEntity +end + +type Camera3D <: GraphicEntity +end + +type WorldCoordinateSystem3D <: GraphicEntity +end + function draw!(g::EpipolarLineGraphic, l::AbstractVector, dim::Tuple{<:Number,<:Number}, p::RecipesBase.AbstractPlot{<:RecipesBase.AbstractBackend}) @@ -46,3 +58,69 @@ function is_inbounds(pt::AbstractVector, dim::Tuple{<:Number,<:Number}) nrow, ncol = dim pt[1] >= -1.5 && pt[1] < ncol+1.5 && pt[2] >= -1.5 && pt[2] <= nrow + 1.5 end + + +function draw!(g::LineSegment3D, 𝐨::AbstractArray, 𝐩::AbstractArray, col::Symbol, p::RecipesBase.AbstractPlot{<:RecipesBase.AbstractBackend}) + x = [𝐨; 𝐩][:,1] + y = [𝐨; 𝐩][:,2] + z = [𝐨; 𝐩][:,3] + Plots.path3d!(x,y,z, w = 2,grid = false, box = :none, legend = false, linecolor = col) +end + +function draw!(g::LineSegment3D, 𝐨::AbstractVector, 𝐩::AbstractVector, col::Symbol, p::RecipesBase.AbstractPlot{<:RecipesBase.AbstractBackend}) + draw!(LineSegment3D(), 𝐨', 𝐩', col, p) +end + +function draw!(g::PlaneSegment3D, 𝐩₁::AbstractArray, 𝐩₂::AbstractArray, 𝐩₃::AbstractArray, 𝐩₄::AbstractArray, col::Symbol, p::RecipesBase.AbstractPlot{<:RecipesBase.AbstractBackend}) + draw!(LineSegment3D(), 𝐩₁, 𝐩₂, col, p) + draw!(LineSegment3D(), 𝐩₂, 𝐩₃, col, p) + draw!(LineSegment3D(), 𝐩₃, 𝐩₄, col, p) + draw!(LineSegment3D(), 𝐩₄, 𝐩₁, col, p) +end + +function draw!(g::WorldCoordinateSystem3D, scale, p::RecipesBase.AbstractPlot{<:RecipesBase.AbstractBackend}) + πžβ‚ = [1, 0, 0] + πžβ‚‚ = [0, 1, 0] + πžβ‚ƒ = [0, 0, 1] + 𝐨 = [0, 0, 0] + + # Draw the world coordinate axes. + draw!(LineSegment3D(), 𝐨, 𝐨 + scale*πžβ‚, :red, p) + draw!(LineSegment3D(), 𝐨, 𝐨 + scale*πžβ‚‚, :green, p) + draw!(LineSegment3D(), 𝐨, 𝐨 + scale*πžβ‚ƒ, :blue, p) +end + +function draw!(g::Camera3D, 𝐊::AbstractArray, 𝐑::AbstractArray, 𝐭::AbstractArray, scale, p::RecipesBase.AbstractPlot{<:RecipesBase.AbstractBackend}) + + # Origin of the world coordinate system. + πžβ‚ = [1, 0, 0] + πžβ‚‚ = [0, 1, 0] + πžβ‚ƒ = [0, 0, 1] + 𝐨 = [0, 0, 0] + + # Initial camera imaging plane. + 𝐩₁ = [-125, 125, -50] + 𝐩₂ = [125, 125, -50] + 𝐩₃ = [125, -125, -50] + 𝐩₄ = [-125, -125, -50] + + # Initial camera center. + 𝐜 = [0.0, 0.0, 0.0] + 𝐜 = 𝐑*𝐜 + 𝐭 + Plots.plot!([𝐜[1]],[𝐜[2]],[𝐜[3]],seriestype = :scatter, ms=1, grid = false, box = :none, legend = false, markercolor=:Red) + + draw!(PlaneSegment3D(), 𝐑*𝐩₁ + 𝐭, 𝐑*𝐩₂ + 𝐭, 𝐑*𝐩₃ + 𝐭, 𝐑*𝐩₄ + 𝐭, :black, p) + + # Connect camera center with corners of plane segment. + draw!(LineSegment3D(), 𝐜, 𝐑*𝐩₁ + 𝐭, :black, p) + draw!(LineSegment3D(), 𝐜, 𝐑*𝐩₂ + 𝐭, :black, p) + draw!(LineSegment3D(), 𝐜, 𝐑*𝐩₃ + 𝐭, :black, p) + draw!(LineSegment3D(), 𝐜, 𝐑*𝐩₄ + 𝐭, :black, p) + + # Draw camera coordinate axes for the first camera. + draw!(LineSegment3D(), 𝐜, (𝐑*scale*πžβ‚ + 𝐭), :red, p) + draw!(LineSegment3D(), 𝐜, (𝐑*scale*πžβ‚‚ + 𝐭), :green, p) + draw!(LineSegment3D(), 𝐜, (𝐑*scale*πžβ‚ƒ + 𝐭), :blue, p) + + +end diff --git a/src/estimate/estimate_twoview.jl b/src/estimate/estimate_twoview.jl index 19603ed..820bc8f 100644 --- a/src/estimate/estimate_twoview.jl +++ b/src/estimate/estimate_twoview.jl @@ -1,5 +1,5 @@ function estimate(entity::FundamentalMatrix, method::DirectLinearTransform, π’Ÿ::Tuple{AbstractArray, Vararg{AbstractArray}}) - β„³, β„³ΚΉ = collect(π’Ÿ) + β„³, β„³ΚΉ = π’Ÿ N = length(β„³) if (N != length(β„³ΚΉ)) throw(ArgumentError("There should be an equal number of points for each view.")) diff --git a/src/noise/ModuleNoise.jl b/src/noise/ModuleNoise.jl new file mode 100644 index 0000000..c6bb26c --- /dev/null +++ b/src/noise/ModuleNoise.jl @@ -0,0 +1,6 @@ +module ModuleNoise +using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleMathAliases, MultipleViewGeometry.ModuleOperators +using StaticArrays +export perturb +include("perturb.jl") +end diff --git a/src/noise/perturb.jl b/src/noise/perturb.jl new file mode 100644 index 0000000..d9aef1e --- /dev/null +++ b/src/noise/perturb.jl @@ -0,0 +1,15 @@ +# Assume homogeneous coordinates +function perturb(noise::GaussianNoise, Οƒ::Real, π’Ÿ::Tuple{AbstractArray, Vararg{AbstractArray}}) + 𝓔 = deepcopy(π’Ÿ) + S = length(𝓔) + for s = 1:S + β„³ = 𝓔[s] + N = length(β„³) + for n = 1:N + 𝐦 = β„³[n] + D = length(𝐦) + 𝐦 .= 𝐦 + push(Οƒ*SVector(randn((D-1,1))...),0.0) + end + end + 𝓔 +end diff --git a/src/triangulation/ModuleTriangulation.jl b/src/triangulation/ModuleTriangulation.jl new file mode 100644 index 0000000..14b0532 --- /dev/null +++ b/src/triangulation/ModuleTriangulation.jl @@ -0,0 +1,9 @@ +module ModuleTriangulation +using MultipleViewGeometry.ModuleMathAliases, MultipleViewGeometry.ModuleDataNormalization +using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleProjection +using MultipleViewGeometry.ModuleConstraints, MultipleViewGeometry.ModuleConstruct +using MultipleViewGeometry.ModuleOperators +using StaticArrays +export triangulate +include("triangulate.jl") +end diff --git a/src/triangulation/triangulate.jl b/src/triangulation/triangulate.jl new file mode 100644 index 0000000..c8dd2ab --- /dev/null +++ b/src/triangulation/triangulate.jl @@ -0,0 +1,56 @@ +function triangulate(method::DirectLinearTransform, 𝐅::AbstractArray, π’Ÿ::Tuple{AbstractArray, Vararg{AbstractArray}}) + β„³, β„³ΚΉ = π’Ÿ + 𝐏₁, 𝐏₂ = construct(ProjectionMatrix(),𝐅) + N = length(β„³) + 𝒴 = Array{Point3DH}(N) + for n = 1:N + 𝐦 = β„³[n] + 𝐦ʹ = β„³ΚΉ[n] + 𝐀 = [ (𝐦[1] * 𝐏₁[3,:] - 𝐏₁[1,:])' ; + (𝐦[2] * 𝐏₁[3,:] - 𝐏₁[2,:])' ; + (𝐦ʹ[1] * 𝐏₂[3,:] - 𝐏₂[1,:])' ; + (𝐦ʹ[2] * 𝐏₂[3,:] - 𝐏₂[2,:])' ] + # 𝐀₁ = vec2antisym(𝐦)*𝐏₁ + # 𝐀₂ = vec2antisym(𝐦ʹ)*𝐏₂ + # @show typeof(𝐀₁) + # @show size(𝐀₁) + # @show size(𝐀₂) + # 𝐀 = vcat(𝐀₁,𝐀₂) + # @show size(𝐀) + # Ξ», f = smallest_eigenpair(Symmetric(Array(𝐀'*𝐀))) + # @show Ξ» + # 𝒴[n] = Point3DH(𝑛(f)) + + U,S,V = svd(Array(𝐀)) + 𝒴[n] = Point3DH(𝑛(V[:,4])) + end + 𝒴 +end + +function triangulate(method::DirectLinearTransform, 𝐏₁::AbstractArray, 𝐏₂::AbstractArray, π’Ÿ::Tuple{AbstractArray, Vararg{AbstractArray}}) + β„³, β„³ΚΉ = π’Ÿ + N = length(β„³) + 𝒴 = Array{Point3DH}(N) + for n = 1:N + 𝐦 = β„³[n] + 𝐦ʹ = β„³ΚΉ[n] + 𝐀 = [ (𝐦[1] * 𝐏₁[3,:] - 𝐏₁[1,:])' ; + (𝐦[2] * 𝐏₁[3,:] - 𝐏₁[2,:])' ; + (𝐦ʹ[1] * 𝐏₂[3,:] - 𝐏₂[1,:])' ; + (𝐦ʹ[2] * 𝐏₂[3,:] - 𝐏₂[2,:])' ] + # 𝐀₁ = vec2antisym(𝐦)*𝐏₁ + # 𝐀₂ = vec2antisym(𝐦ʹ)*𝐏₂ + # @show typeof(𝐀₁) + # @show size(𝐀₁) + # @show size(𝐀₂) + # 𝐀 = vcat(𝐀₁,𝐀₂) + # @show size(𝐀) + # Ξ», f = smallest_eigenpair(Symmetric(Array(𝐀'*𝐀))) + # @show Ξ» + # 𝒴[n] = Point3DH(𝑛(f)) + + U,S,V = svd(Array(𝐀)) + 𝒴[n] = Point3DH(𝑛(V[:,4])) + end + 𝒴 +end diff --git a/src/types/ModuleTypes.jl b/src/types/ModuleTypes.jl index e4ca9d7..1bf3158 100644 --- a/src/types/ModuleTypes.jl +++ b/src/types/ModuleTypes.jl @@ -2,14 +2,14 @@ module ModuleTypes using StaticArrays export HomogeneousPoint, ProjectiveEntity, FundamentalMatrix, ProjectionMatrix -export HomogeneousCoordinates +export HomogeneousCoordinates, EssentialMatrix export CameraModel, Pinhole, CanonicalLens -export EstimationAlgorithm, DirectLinearTransform, Taubin,FundamentalNumericalScheme +export EstimationAlgorithm, DirectLinearTransform, Taubin, FundamentalNumericalScheme export CostFunction, ApproximateMaximumLikelihood, AML export CoordinateSystemTransformation, CanonicalToHartley, HartleyToCanonical export CovarianceMatrices export Point2DH, Point3DH export HessianApproximation, CanonicalApproximation, CovarianceEstimationScheme - +export NoiseModel, GaussianNoise include("types.jl") end diff --git a/src/types/types.jl b/src/types/types.jl index 816bca3..4506044 100644 --- a/src/types/types.jl +++ b/src/types/types.jl @@ -47,9 +47,15 @@ abstract type CoordinateSystemTransformation end abstract type CovarianceEstimationScheme end +abstract type NoiseModel end + + type FundamentalMatrix <: ProjectiveEntity end +type EssentialMatrix <: ProjectiveEntity +end + type ProjectionMatrix <: ProjectiveEntity end @@ -77,8 +83,6 @@ end type HessianApproximation <: CovarianceEstimationScheme end - - type Taubin <: EstimationAlgorithm end @@ -94,3 +98,8 @@ end type CovarianceMatrices end + + +type GaussianNoise <: NoiseModel + +end diff --git a/test/construct_essentialmatrix_tests.jl b/test/construct_essentialmatrix_tests.jl new file mode 100644 index 0000000..a89d03a --- /dev/null +++ b/test/construct_essentialmatrix_tests.jl @@ -0,0 +1,28 @@ +using MultipleViewGeometry, MultipleViewGeometry.ModuleRotation, Base.Test + +# Intrinsic and extrinsic parameters for the first camera. +πŠβ‚ = zeros(3,3) +πŠβ‚[1,1] = 10 +πŠβ‚[2,2] = 10 +πŠβ‚[3,3] = 1 +𝐑₁ = rotxyz(deg2rad(10), deg2rad(15), deg2rad(45)) +𝐭₁ = [-250.0, 0.0, 2500.0] + +# Intrinsic and extrinsic parameters for the second camera. +πŠβ‚‚ = zeros(3,3) +πŠβ‚‚[1,1] = 5 +πŠβ‚‚[2,2] = 5 +πŠβ‚‚[3,3] = 1 +𝐑₂ = rotxyz(deg2rad(10), deg2rad(15), deg2rad(45)) +𝐭₂ = [250.0, 0.0, 2500.0] + +𝐅 = construct(FundamentalMatrix(),πŠβ‚,𝐑₁,𝐭₁,πŠβ‚‚,𝐑₂,𝐭₂) +𝐄 = construct(EssentialMatrix(),𝐅, πŠβ‚, πŠβ‚‚) + +# Result 9.17 of R. Hartley and A. Zisserman, β€œTwo-View Geometry,” Multiple View Geometry in Computer Vision +# A 3 by 3 matrix is an essential matrix if and only if two of its singular values +# are equal, and the third is zero. +U, S , V = svd(𝐄) + +@test isapprox.(S[1], S[2]; atol = 1e-14) +@test isapprox.(S[3], 0.0; atol = 1e-10) diff --git a/test/perturb_tests.jl b/test/perturb_tests.jl new file mode 100644 index 0000000..9a58491 --- /dev/null +++ b/test/perturb_tests.jl @@ -0,0 +1,23 @@ +using MultipleViewGeometry, Base.Test +using MultipleViewGeometry.ModuleCostFunction +using MultipleViewGeometry.ModuleTypes +using MultipleViewGeometry.ModuleConstraints +using MultipleViewGeometry.ModuleConstruct +using MultipleViewGeometry.ModuleNoise +using BenchmarkTools, Compat +using StaticArrays + +# Fix random seed. +srand(1234) + +𝒳 = [Point3DH(x,y,z,1.0) + for x=-1:1:10 for y=-1:1:10 for z=-1:1:10] +π’Ÿ = perturb(GaussianNoise(), 1.0, tuple(𝒳) ) +𝒳ʹ = π’Ÿ[1] + +N = length(𝒳ʹ) +for n = 1:N + @test !isapprox(sum(abs.(𝒳[1][1:3]-𝒳ʹ[1][1:3])/4), 0.0; atol = 1e-12) + # No noise should have been added to the last coordinate. + @test isapprox(sum(abs.(𝒳[1][4]-𝒳ʹ[1][4])), 0.0; atol = 1e-12) +end diff --git a/test/runtests.jl b/test/runtests.jl index e952379..34da0b8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,8 +6,12 @@ using Base.Test @testset "Data Normalization Tests" begin include("data_normalization_tests.jl") end @testset "Transform Tests" begin include("transform_tests.jl") end @testset "Fundamental Matrix Construction Tests" begin include("construct_fundamentalmatrix_tests.jl") end +@testset "Essential Matrix Construction Tests" begin include("construct_essentialmatrix_tests.jl") end @testset "Estimate Two View Tests" begin include("estimate_twoview_tests.jl") end @testset "Projection Matrix Construction Tests" begin include("construct_projectionmatrix_tests.jl") end @testset "Rotations Tests" begin include("rotations_tests.jl") end @testset "Cost Function Tests" begin include("cost_functions_tests.jl") end @testset "Fundamental Matrix Covariance Test" begin include("fundamental_matrix_covariance_test.jl") end +@testset "Satisfy Epipolar Constraints Test" begin include("satisfy_epipolar_constraints_tests.jl") end +@testset "Triangulate Test" begin include("triangulate_tests.jl") end +@testset "Noise Test" begin include("perturb_tests.jl") end diff --git a/test/satisfy_epipolar_constraints_tests.jl b/test/satisfy_epipolar_constraints_tests.jl new file mode 100644 index 0000000..347588d --- /dev/null +++ b/test/satisfy_epipolar_constraints_tests.jl @@ -0,0 +1,86 @@ +using MultipleViewGeometry, Base.Test +using MultipleViewGeometry.ModuleCostFunction +using MultipleViewGeometry.ModuleTypes +using MultipleViewGeometry.ModuleConstraints +using BenchmarkTools, Compat +using StaticArrays + +# Fix random seed. +srand(1234) + +# Test for cost functions. + +# Test cost function on Fundamental matrix estimation. + +𝒳 = [Point3DH(x,y,z,1.0) + for x=-1:0.5:10 for y=-1:0.5:10 for z=2:-0.1:1] + +# Intrinsic and extrinsic parameters of camera one. +πŠβ‚ = @SMatrix eye(3) +𝐑₁ = @SMatrix eye(3) +𝐭₁ = @SVector [0.0, 0.0, -10] + +# Intrinsic and extrinsic parameters of camera two. +πŠβ‚‚ = @SMatrix eye(3) +𝐑₂ = @SMatrix eye(3) #SMatrix{3,3,Float64,9}(rotxyz(pi/10,pi/10,pi/10)) +𝐭₂ = @SVector [10.0, 10.0, -10.0] + +# Camera projection matrices. +𝐏₁ = construct(ProjectionMatrix(),πŠβ‚,𝐑₁,𝐭₁) +𝐏₂ = construct(ProjectionMatrix(),πŠβ‚‚,𝐑₂,𝐭₂) + +# Set of corresponding points. +β„³ = project(Pinhole(),𝐏₁,𝒳) +β„³ΚΉ = project(Pinhole(),𝐏₂,𝒳) + +𝐅 = construct(FundamentalMatrix(),πŠβ‚,𝐑₁,𝐭₁,πŠβ‚‚,𝐑₂,𝐭₂) + +# Verify that the algorithm returns the correct answer when the +# constraint is already satisfied. +π’ͺ ,π’ͺΚΉ = satisfy(FundamentalMatrix(), EpipolarConstraint(), 𝐅, (β„³, β„³ΚΉ)) + +# Verify that the original corresponding points satisfy the epipolar constraint. +N = length(β„³) +for n = 1:N + 𝐦 = β„³[n] + 𝐦ʹ = β„³ΚΉ[n] + @test isapprox(𝐦'*𝐅*𝐦ʹ, 0.0; atol = 1e-14) +end + +# Verify that the 'corrected' points satisfy the epipolar constraint. +N = length(β„³) +for n = 1:N + 𝐦 = π’ͺ[n] + 𝐦ʹ = π’ͺΚΉ[n] + @test isapprox(𝐦'*𝐅*𝐦ʹ, 0.0; atol = 1e-14) +end + +# Perturb the original corresponding points slightly so that they no-longer +# satisfy the epipolar constraint. +N = length(β„³) +Οƒ = 1e-7 +for n = 1:N + β„³[n] = β„³[n] + SVector{3}(Οƒ * vcat(rand(2,1),0)) + β„³ΚΉ[n] = β„³ΚΉ[n] + SVector{3}(Οƒ * vcat(rand(2,1),0)) + 𝐦 = β„³[n] + 𝐦ʹ = β„³ΚΉ[n] + @test abs(𝐦'*𝐅*𝐦ʹ) > 1e-12 +end + + +# Verify that the algorithm returns the correct answer when applied +# to sets of correspondences that do not satisfy the epipolar constraint. +π’ͺ ,π’ͺΚΉ = satisfy(FundamentalMatrix(), EpipolarConstraint(), 𝐅, (β„³, β„³ΚΉ)) + +# Verify that the 'corrected' points satisfy the epipolar constraint. +N = length(β„³) +for n = 1:N + 𝐦 = π’ͺ[n] + 𝐦ʹ = π’ͺΚΉ[n] + @test isapprox(𝐦'*𝐅*𝐦ʹ, 0.0; atol = 1e-14) +end + +# Οƒ = 1e-7 +# SVector{3}(Οƒ * vcat(rand(2,1),0)) +# +# β„³[1] - π’ͺ[1] diff --git a/test/triangulate_tests.jl b/test/triangulate_tests.jl new file mode 100644 index 0000000..9a27898 --- /dev/null +++ b/test/triangulate_tests.jl @@ -0,0 +1,62 @@ +using MultipleViewGeometry, Base.Test +using MultipleViewGeometry.ModuleCostFunction +using MultipleViewGeometry.ModuleTypes +using MultipleViewGeometry.ModuleConstraints +using MultipleViewGeometry.ModuleConstruct +using BenchmarkTools, Compat +using StaticArrays + +# Fix random seed. +srand(1234) + +𝒳 = [Point3DH(x,y,z,1.0) + for x=-1:0.5:10 for y=-1:0.5:10 for z=2:-0.1:1] + +# Intrinsic and extrinsic parameters of camera one. +πŠβ‚ = @SMatrix eye(3) +𝐑₁ = @SMatrix eye(3) +𝐭₁ = @SVector [0.0, 0.0, -10] + +# Intrinsic and extrinsic parameters of camera two. +πŠβ‚‚ = @SMatrix eye(3) +𝐑₂ = @SMatrix eye(3) #SMatrix{3,3,Float64,9}(rotxyz(pi/10,pi/10,pi/10)) +𝐭₂ = @SVector [10.0, 10.0, -10.0] + +# Camera projection matrices. +𝐏₁ = construct(ProjectionMatrix(),πŠβ‚,𝐑₁,𝐭₁) +𝐏₂ = construct(ProjectionMatrix(),πŠβ‚‚,𝐑₂,𝐭₂) + +# Set of corresponding points. +β„³ = project(Pinhole(),𝐏₁,𝒳) +β„³ΚΉ = project(Pinhole(),𝐏₂,𝒳) + +𝒴 = triangulate(DirectLinearTransform(),𝐏₁,𝐏₂,(β„³,β„³ΚΉ)) + +# Triangulating with the same projection matrices that were used to construct +# (β„³,β„³ΚΉ) should yield the same 3D points as the original 𝒳. +N = length(𝒴) +for n = 1:N + @test isapprox(sum(abs.(𝒳[n]-𝒴[n])/4), 0.0; atol = 1e-12) +end + + +𝐅 = construct(FundamentalMatrix(),πŠβ‚,𝐑₁,𝐭₁,πŠβ‚‚,𝐑₂,𝐭₂) + +# To triangulate the corresponding points using the Fundamental matrix, we first +# have to factorise the Fundamental matrix into a pair of Camera matrices. Due +# to projective ambiguity, the camera matrices are not unique, and so the +# triangulated 3D points will most probably not match the original 3D points. +# However, when working with noiseless data, the projections of the triangulated +# points should satisfy the epipolar constraint. We can use this fact to +# validate that the triangulation is correctly implemented. +𝒴 = triangulate(DirectLinearTransform(),𝐅,(β„³,β„³ΚΉ)) + +𝐐₁, 𝐐₂ = construct(ProjectionMatrix(),𝐅) +π’ͺ = project(Pinhole(),𝐐₁,𝒴) +π’ͺΚΉ= project(Pinhole(),𝐐₂,𝒴) +N = length(π’ͺ) +for n = 1:N + 𝐦 = π’ͺ[n] + 𝐦ʹ = π’ͺΚΉ[n] + @test isapprox(𝐦'*𝐅*𝐦ʹ, 0.0; atol = 1e-14) +end