From cfa872fcb5943d8b8e0a822fa0a3d0afc2d9c9a0 Mon Sep 17 00:00:00 2001 From: Jonas Ries Date: Fri, 4 Dec 2020 12:34:33 +0100 Subject: [PATCH 1/4] Gauss for BP --- other/BP_Gauss_analysis.m | 74 ++++++++++++++++++++++++++++++++++++++ settings/cameras.mat | Bin 18664 -> 18689 bytes 2 files changed, 74 insertions(+) create mode 100644 other/BP_Gauss_analysis.m diff --git a/other/BP_Gauss_analysis.m b/other/BP_Gauss_analysis.m new file mode 100644 index 00000000..08f6db30 --- /dev/null +++ b/other/BP_Gauss_analysis.m @@ -0,0 +1,74 @@ + +pixs=100; +dz=10; + +%% BEADS: calibrate sx +figure(88); + +sx=g.locData.loc.PSFxnm/pixs; +sy=g.locData.loc.PSFynm/pixs; +frame=g.locData.loc.frame; +z=frame*dz; +plot(frame,sx,'.',frame,sy,'.') + +ds=sx.^2-sy.^2; +figure(89) +plot(frame,ds,'.') + +range=[45,95]; + +inr=frame>range(1)&framemaxrange(2); +g.locData.loc.znm(outofrange)=outz; + +%% data: average x,y +trafo=g.locData.files.file.transformation; +inref=trafo.getRef(g.locData.loc.xnm,g.locData.loc.ynm); +cr=[g.locData.loc.xnm(inref),g.locData.loc.ynm(inref)]; +ct=[g.locData.loc.xnm(~inref),g.locData.loc.ynm(~inref)]; +ctt=trafo.transformToReference(2,ct/pixs)*pixs; +lr.x=g.locData.loc.xnm;lr.y=g.locData.loc.ynm;lr.frame=g.locData.loc.frame; +ft=g.locData.loc.frame(~inref); +locpt=g.locData.loc.locprecnm(~inref); +lt.x=ctt(:,1);lt.y=ctt(:,2);lt.frame=ft; +[iA,iB]=matchlocsall(lr,lt,0,0,500); + +normw=1./g.locData.loc.locprecnm(iA).^2+1./locpt(iB).^2; +xav=(lr.x(iA)./g.locData.loc.locprecnm(iA).^2+lt.x(iB)./locpt(iB).^2)./normw; +yav=(lr.y(iA)./g.locData.loc.locprecnm(iA).^2+lt.y(iB)./locpt(iB).^2)./normw; + +figure(91) +plot(lr.y(iA),lt.y(iB),'.') + + +% plot(g.locData.loc.xnm(inref),g.locData.loc.ynm(inref),'.',g.locData.loc.xnm(~inref),g.locData.loc.ynm(~inref),'.') +figure(92) +plot(g.locData.loc.xnm(inref),g.locData.loc.ynm(inref),'.',cn(:,1),cn(:,2),'.',xav,yav,'.') + + +ind=true(size(g.locData.loc.xnm)); +ind(iA)=false; +g.locData.removelocs(ind) +g.locData.loc.xnm=xav; +g.locData.loc.ynm=yav; + +g.locData.regroup; diff --git a/settings/cameras.mat b/settings/cameras.mat index 9dee023961b3b7ceab8d83f3980674d9a8c884f5..60f284e153d707ad41139f77029724646a43ff7d 100644 GIT binary patch delta 16871 zcmXtfcTiK&_BBPMh;#&{iwJ^rr9)Ii6a+*>q)QV4DWTU~1*8}0od_r>2uSY? zBcmiM$1fu#BbA^pnDG8;g(rv-(9$3y^UdDLcxcJ0 zp5r4egfRq~KhS-)snK2)(6$dCeW@dJP#0}p_iZ~!YTI9hn=Fw=u>OCQ& z2=L3^M!tQ~W6y`GTF;i~t*wy;DD?C6*cFR76lIKcfp-8usw{EX{jp!?H*ra~2$dcD zUBZ{{@^zu>>QApo`2j((fRRjOyLKn+j4wTVaqbfVd%4s_n5s1O>Ha?b(zIFnp~eb7 z%IdkQttx{-A7>_=Vtf!cf0N(pie3LaNi;J5Mk=@wX_8$FM3AyX4j8{Zsc8Y z?wI})&Q0a6SZS$edU0`(M0b(bEjDeriw@)cAog$qujHkofJ!8}^q@`^C|)5Sa#1%@ zdE2h&YIdL4LkE`QV7%GflZZ|4LZP)0%)Z>*kM_ZQzdovf-Qf+HN0NhI72A;Y` z{;S54Z%2k8L{}iUNT{UfH*xC3cgo$%yb9bXs>Zv%(|A$k5a_iw|KmH#3+(Xwp)PM( zLw>W(gl>xKOCnD3xlAniI)!gfd{X=i2KQ1Rx1#^|kZYDJu6T#D8(ooWx?n-QI$3mg zhthihRbz>Ih%GQx<9qXG6s?|xb-^rTbwlO$iLMI)1^|fVzD%Ff)cR63w>=aqJ5V(3 z<}m%>!g&-fv-c-Bpm_SX?IwXD=_VhGXjuCDFIN15S(Ba3nc_yNM%#*6y9KXF1}Nu9 z&w!g>7hf(|w(&pd^{9oUud4erg3gFP0_S%_P6IrOHG2VaGBOfr#~7(ME5+vBA5U%_ zqr>_s4G=HV|DiAB?9gBO(KLu;%v8ss=es&Y!yOGz0Gu(>hNK zSZ_?;OV)ZCHDYw82FK#Y_Z7pujWjtwvy(kaTEHAze8=9WDR1bJcnXL;MU_9RwaHYw zBVaMmilkq%|BU9(19P+<)c*_PC|jZCP3~f~1=J}s3TZBw!|lBh70wYaP^7n4AzwRW zxu{QX9y1JJb6wOeUy>B#!l$7LkNqD(E=~ouP9r$gbf}fWYi%Am@F>Nh+ZY}*$R$v9 z5uRu1ak0z~MM@-v-)&?%BGJlAR0ZX;{^A_^Iu(HtC?{dN;$N)Tp}`{%J<+|{_5p?^ z-~;&~MV^_i#2~ak{O@Sl{66_dbp!8 z|8r{JF=xdXzC6AA?NJ?r)R%<@7w74VoDEERhwSVQHiI&oI4a{E+%E5Gfm!N{*=dR0S)?V zo(UP{-MTF^Y&#rfn?m}9$A31O|B^-I(AfVgrVR6-Z+>D**i5~8dIn);2xV?qG3|6(?^?LXYq zol`1C*)}Hd0Q@gOyTq2qDM^Qp1CYHnDiPO_M`K!9EPR5`#T*fiwgni!(&Vu=2)p~N zsiiiBKuQa|c6y2a#FDqJJZSzah0$|GFVE7gzZ9^c)c+8IZ-L5}(WMvg_(olg$;Zk1 zuNvj4;a<6gv;QI|{JJ0d<$gRdJmg@eJ?Kb1$E2QHX%jWHiLF05uuS+vl5QWH+kb;m4m!Zo>;I#Zwv&FKfXZ0+IbHNI zi1@=K>(|L9%vnE8+;D`54hnp2_&Zrb*SlRa0xuZ`jUJ0KXom395H?pAnkC8~@UdH# zAQofJM{lz$XahZYCg-Yv9JH@y(`YMoex7(rL7aBgK`)*p>nAmyRAYXt#`oOs@cFiW zuwK&W$VYKPHsL1eZSjbgR_sdtVHKE6eCL|rg}!(xwk3mh%jJIuR}f_f@epJCBYM@!Ogta2S>ty0F*tbRBFjRd|SupK8)Uiv2I3Jl)dI|@_a0|^F~OjiX?#@T~% zoggHo-;92fAF!_A3-XK$;h|Q({9AKuNQdg`qA8*+F2w0LE`%h>a6$Xu5=w9TPuuS{ ze@pk0OO$&kPSxFE^}+Uri~4A-g=pdKUz^^{w29ZBjN)VkF3zW&^Tv|9rZPiCO)K+# z&VrDB5O1}gv6&)~S2++g0)IYe52!!P`lxtn1P0Uxj7o}lo!PrM1v$#XS4{R?{xq!p z?VX%jyd`qF8uoBb*(OIx;;~kHpjVXQxhf52BrQzgXy&;Gw5<)+D|;??Fvo0`tx z)hnU%xQ{pel_c^$M|z^YzK$KfJAgW#35HC9xvgT93zl|68pk%9r9C|f8C9UxJ4^S) z_dLzuTvZdJL8muV_o~nffaSC7>1O-ccK|duDLD?kdNt~GkWe>!&9I;I#)fbHwvd)# z`AIB(!{RN4ncoW$+dD|}{TkY$Hzy^A za6Yc)IAaYUq`nUxmHfoPh!mu*iGV-9Yu^(@3VT#2+mE*~@1N0-$o6mmpEHo3<+UmQ1i`V)V27|2xCxFI9t4>cD6W= zZ7W2*_+|sm8l&@GJIO=C-88@{+C-n>>o99+$x^|5074-_Xhr*lt`$5e)fSn5V+ihg=!aE-`7DligGJIz0TvzB&?i# z@b-@m9lQE_@8x+4aM|nbv$kr8+M0EPhSaNRF`i~?LmB0XOMqzWy&f-*dSn!Ign=2I zPOl#Q%Ec;JF-#5iv@DXZHjkt?@Z0aHNF>I)LTO2qE-jap40Cc4tsZy3h#w`|BgwI9 zM0Zp|g{uh;Y~n{n;*zM3G`0*M%UK>xtZ3S)Ik+;VRsIYLOO*>A9%1NwYUye*oCP1o zM}F59Tu4%G0bmQk`E_mubsUHKo%nVQz5q_?mJ_fMkuD+rYUatdgvGPxKYad`qYqr! zEoS<4W@}3=urlmvDmjc%an?fUv_`(rdrbku1MBG;2UiElCs2 zvi!>B)x%v8Va)1f|r7RjBS$hmg~qEL$3R{s;;+ zhg|)KLOJaLt|t`(?6PG$K>r0cEbAEVZP#w-v0idb!WgL36R@6+Hm70e!YLjXcv@Oj zLQl8>F&FWd+xuRYmG_#hr4DUCnuynz3)(H&X#G)#%6>FAHVw~l-Mzq~maA8)O(#;u zoC6!-YW?`lpu#^$yC8{e4E#6pE|~LGcPAr#0)qE6!I;|58d>W}E2R#rE4+}Qz{VmgB9r4F%Z&(&_JZDfSITN@CbvcGxPApbQvCvQgkses+{}wPp{)f}O z8*4wQzsVZ9qbssS%TSW&^7|@>|GneaziK<`bsSbaQ70B(2XqS!V!)Uu&lEoXsS?LB~0q06YwO=eW(@EEE)~(||U=XJ(Mh(=lnYZ># z3_5twt>QDHy61DTnf&`HjxtK~J9r2%q5D4!sgqMOcoxaAqdMz^NRrN7 z zl$rpc_unRv$y;%is64+D!(>_-A4P)FcGVjF*(3fz!Q-Dgr(*TrW*0|-i@(eB!yL`I zX3suEtq)xr;%TO>aBN0&LLcAaUoKzIV{wfT?)Uf}EL_IbnvHjs{DPuP{%*Fi7qVL4 z9)`I``-IN$klt#8f6{Qk)FxQ^IOE6pU2kbb}femagz>lUw@MlawHv>Z9@5YTul?``h%tc8x}*k`bx_KEqq9oIb;~M8>IP zq#@YY>o8@t=&zs5`yS|(M!od7z|L>GhcD5;&a6bvE!idqT7MGvHBRuJ`Y;3tE)*GFldnFT8zC1SkS2H*dY2GZ$2 zq_HM9(D+ze@+%L9acrv02ADCGn0D`#YwU>YiiH|{5n?xG$SWc$5@qehPgzQ66W7}0 z?tXb;(SvR&!XB?ghtuSZZr%Jopzzzd$Up3)h?Aq{y6sJys3E@-V&O&UHRMP&8Ei?F z3Yc4t*B|`*cFtfdPVCb=P(}W~wq$hrqu|ju;eg;_yM@$~&qq+49>7Lu=)9Jdx=Mb8 z+zeN_`)XsH^z=w@$5mk|zNZw2!o zDxfei^-{0U(=qdnQF*LL2Qb4?fV(Ue)}qf~YJ44?E)mUa7Z`ePU19J*_Rd{D-!L3@ z?JnJ3=>u0QJ)`@kK-;*d1xbR?tY8X?TP>On7#t+&oX#5YdT`=Yao>f3`6R}Gqm|si z>P`Ic5`-4e{54-5d3yoo0jbjMhPa&^PO6s7*rz2z`%|IH5v*^InjCrILgkY8Csjz; zc|I!JzeCNWH6q8yb&_tB*v?_#smf6S^}jlzbT5gu{{@i!8{^sy^Ek&J94INNmrv?X z2DmHPzbVq?cnw)L?p(c~Ef)|gHo(E|6rNZ}H3m1g#<-xZUb~jamo5;DP6}VN1WsP_ z@|qTycVNMupLVy5N6gZffJhcZzqWrWAQr>L@Ridh9 z%c`=%uGR^aU1O8JZ#>|1-hyob?C~vWZz#!B-WyOJrrWRx(fq62mj9$%rP%QS7625B zIsHa*a`Yb@hhAa+T71TY8+u$O76OPYHq7zE7{OCZmX}5bov*}74=C)oph?4-@9HB_ z27_Q;3E;8r6ovbl137YyCnB@-6$40w;nVe?$8Lq|oHoneGH%E#xn>OSId9L=N2P~qgz86p zVhs8kN5^QMK&r2Iyt-GhD4*u!yj-)?$D0S+6*{fS%{2=n-urM*F%9#{I+fv_cV*U2 zqjLcB&V_E6<>N)fnt2r`VDNIY*94{U7`fa4O|>c;6R>k?n(Vc7Z|0>vuWzd;C$~TAIjx9R5kc4t`JYPx)jZ18hNpSax+$lw z!C7o5@a}~xIQPLt;fmL2P_@>D%)i&u6pkl?n?;&kDLNQfJnt98kMj+EtmNHyjwFT- z9@|9=8u0p2_eHqQ+qU24@st+~jX57#1K-Zhe<&SC*W=dDO>ds-Jcj! zIC~&9R{B`Z4>=n1i-%>#2czoUYg{_Jx+e%C&OO|$eo=T})bq1lxr6gSy)cWqsAiIk4B#LcV1(M8>Js2OhYa8>wMAj-MYj zl=xq5K`i;Zts=BOo5AHJH!?VI&DJ^Kt)MAGoR2y2NNUdT5H@hpzf9KvS=fXrBFY{N zS*CTW2rS#v&F@bWp#VTKHtK51sdzQeXAt5mPtN&{X-~~<Ch#ROj-q9Sq?!M@8&Ypm_#aQvT#><=b>_svDj=z)jn^(-T3P40(}+;I z_vcU{Va4BYBm#j@LpU&#Nr`STeN0(;9;YqPsP_y@s}JHJ7kRczN9STT_N5kE*+QFm z>InL}VLpOnnks&M8(ea;{)Gm}*l9@FdHLGv`I(p!aJ^YYMURd5-lT}4=08O8ot1aJ zi}yrTUO%(ovCfY#idzxPMP3fhXOCM&2FLeBHg=5qHu7H(KMR3286Q>wIsUUhSW;C@5WxebESbRS)ra+ zs^?zQr!Nfc`dvhO+^EM-;9fa%RwZtLQ7mO`^dug$)}^%><$A1vt(O2?QBdBMa~m@t zX1gVFTgm&Ee)lV#JuU`seKpmQ^-*>1nQ!VmAkPC02_p0!QEM9gO|avc&eNSKK9>>_ z5H!3}cG!sVo%Bbg;U;UKE|A=xSOX2$d(qr2)t^Eq4_b6dsSU4iP@F$b+7fQ$=TnaH zGM{t2Dj|O1oMV9HPAMHBFxx4vxS4Qc6U5p-x*jN9B*qlns^^wBkC}?=xgU_HmufzG zwBi#U5Eg6PL;eT`*14x*z_T(FB{K+w?u&%Dnk;4QEt-}bk~Ic|3f^V##^eHBv#tTL zV!rp3t#~ujcLY`od|DVr27gvyC%^x5h|M{0W1m)UK|R<@&aY=(@ghtYzoGAV$HLKX zLMKTJ9Snp&AujDG2;CIM8nVw`_G@vOWo3C^t*#B&JX?{6 zR%K3RtNO`$Jj$S+&u9>Xjnl)bgA9jKI;Xu(Z5dYAU(CARiPs5{<*50^dmE`<6YsL% zkLm_bo^DXNd|GT=*v0L9_Kot;Bx^)4dhuNvr=^PlJRp~1q?4Pfi|tZ(J{AN-9F;Ws za|HwrcOI%OFOT%s9}>smq^H2gc-JFRxUXdrh81GuR)PEyn|;@|QpdAtnKW8x^mr-0 z;hOWV*HP%o0X}S;oy759qv)N!>x>kWuE>zfXtASs(xLzxQ2$E%MGOnIA+Za8vr6D* zq=p^RHbsu5_cq5-kdK*U8oTF?wdUg%Mk^m z!>&1l^m>Z|QFG7So5f-Or2ZZ|vJwd`Kt}d8O z8FVP_%<*{kmo>4rRJRPH@|5_pb$8H%*o69h8lbtnne5v~QAS7n5Dd=)M}*A70u0sa zHsh^~%f|h663R`=$NhILN1Y0TPD@WRe=`8ZCo`nqleIZ0TYquaw1D4}#cHmzJ}re& zXjf?6=2qQ#bLoU@4P+KXJe@#0SAtegFQ@ui6L@cZ0!2~4!ukwXZdf}I9nc3y=mkD) zraveSMH3*AA#8NbcMtt2?Rl)fXx6~cEXW4KF)=I?rl!SsTpZ((DfXcL$ws>ufd5*r zRE~8z1S)=?nER{<6eB> ziPUfHNXk?eE=Aq1r$VQB-qRd2L2`EphmRP-YCXGpqAVA&N4Ykj*?1YU?{`EI0!M%M z_anT`$E@8!jA1{&^(gA^142f|eUWp8RuMGRiLM{dc)q5|M0Q;3)FsyayRQGBxu~k( z`$oc%(axbnPd5jN{S^w@ej3m^4IU*zz{bmmAERp1<) zHG7Yk-0DfvQhYjhVzyTn40pO%sJ$2gH{6quKWlNA1B(v*rmVtk0dL}BcD{xoL*fe7 zMr1aHd(Ko0l+j@E>)59i27S!S+)28&tbNWer4TNk0e;sEaD|H*B&whC39U((tg)mZ zTf|ylU<(m?78V9CKX`rNc|Le5#kKKGV+;3v81rUWa4nT9iZR}LSOFLf403my+aDJB z@bQn&JRY6bu072M>`)`jJ-JT58qLR(**0x0uM4gzzAoc=>pZJgJpPI6iO9^5AI>@w z0`e5qy@fEudxHBiN9?~SJ&2ATVHVV_FLT%Us;%8;>_IRMb^5rI{ZQ}rU@>p3Bt9=OKVYt@BnB-{Cd?OA#O0?>;Yq78ke7h10G$QA8_ekeeq}`&E;=E|2 zQc&v0pa+Xb!7BY(;2L@}lef2^b%^YrZpArNRt~o`bIBFe%;PAy!HWWLvn zvEyx>qgMbOMyx(l!#F1 zbAYCN7U^9LHh1Oxpv|k`#PvQlJzkOh$9)Ba4B)z^zo%clss@rbJtf?UPCEk+d0*ZT zq(kHM5dGc68FF(Pi$5fInLt1{UFEiqKjh?gk*NOL%rpLb*o&BlhP>dfrp)x7x(z0~ z^Xql-3!Nl`kpwuU0f-G9EMHTK0>j zm8di{QvvI&x<@jwUk8nQv8?rf z^&6bWtapwjWjpO2`90aY9~MwLb+2FT5*UT9OnuicficwO(=$;axEYWmUvz5ugf=r2 zf1{&Cf5tdfFG45d3s*8{EOAA>suJ)ARJgD^nqG;Sc zqs)24(@|`gl%8g`vE3K z681y&@@nRW@$qFK0#n`Vx%&O}&4EKdlM6M@%Y;Zk?gjP#@Z{>mD$UvIVJ<)W#p6ti zX>~&^r~>SqL;v{KSnu7KNGeiGhOy-2Rf*nzYu$2(qh}6;D-fQf?Hp5e?UCyFV!3bY zx!c5Cm0lsg+*M90)IS80&su@1dvT*-hBLW-tUG6EeoZZ)Bzx- zPM}+)Hyt-6?>G^?nFFo7v(}Tsk^Y@YGvPX*xsozx7s@XxgEzejQC~@UY2Y4wrhlBN z9?a5oVL5&*f}>>6;H!tjgF*v>tPt^A)Hnv(YY848&ztC9Vp0XnaS5{uQ3p7nts%xE z_zHdAXYGqqoCcTu0u$N}IqSA20F==O@A{9%1~jf%+s}{awpeaGe!m*!NmQliKR4Uu zM(lFcQ)lLAr4GJ(KcgMk7aD)>250vf%U29Hqx`_;N`y>Q;gT8QL$LK`?Hq9r-dRUa zXzUil%@2R!v(DXo^f)n%ec=hnnY z#)IDvufgm-#z$L9P3tEbf!bEhj7#em2=%uI#w6dA*C%#wdY5TwN1tpnzrs+^v8_D) z%%oBDt;GIlIsW3}VUtP-WD?U9SJ}aLcBthL7_`1Fpflqb70a+FZmrQ)697N?; zjX57<(jFYkaI?m^NbVv2&E6CDz$m1C_M*0DHp}E$UjN3;W!ew@`y!X0_tLO}mrC(h z^hyWe$6(e^v%9Ef0T&p$yv!`&t~sgD0i`g=rYV7RuE7|&Rw@*y1ELp&?lb*{G^eX4 zH)9NvLFh(Hza(IU_e=6QRMMK}gFi$;gu%ji)@b}#6c<#~!al|o`vF3d|35wU-y;n~ zP8q!WX$BLuu3@|zSxwTE_n*#k_w+|kC|tUnEdFBxOb`A#i#`>0r)+ZmXg$;ARmY;! z%l1?91ig*6MDCeh6Br8$er2wjwRJkQTr9o|q&9N>#?I{J8{&%f7ffN=eSBiuyHo|4 z+h{k$O&4&y+k@=Ho1wvkP%G}~4HYTRYE|oFx_bVYvKer0fq}ed_4Q#+)~&%RekT6i z^F#g)KK{W7!xo+7QlnX^EE#4Gw>6>2R^i9*JLKtM()%P43ekjtvl@Df*#~hFUjTYw z0nGj7?=~qjuVA;lGmUjs4MNGk5?Q@6c@lD~BLcTL$0*pX81o|&K~CG|%u8ka)*buV z!x+K*e8(m-x|O38+Cq;W=x-XE44U>>wNS)OdmL9y?ykG{rtoUo8beEZ1FGSdRd&x6&y5)by>)pmd^S)^&Q{?u2GX%4DTM^kr{k@jUaj>x@m8kkn( zVY|UC?inU{qnnu>gDt}+Pb}?(i#!4D?oOR##8IC-g#1|1(+HxVE-QXc{Hh%^dE-QA zGOEf6rDMi7ucF4N?;Thzeg|}D9C)+6>eBVkF&#F>Hkoh#s;OtX_p!+gsRw8{AxL|S zMtilnn+iD`Pq6Ri--rjW ze?52?1Z%gV?{}9wTuqQ>n1NugT?3TU*il-#?6!Fq@;m&GJ_5n9_UKHVw}hO zw?}_&yc_u2$XCFkN+%7z>A*w(kp_sWLG=Si4VaMAi#wphu{)#8LdZ z3_Pk+H?yQ?_^irysp(sgM9oVRN@DplJ5{8TuFSX|Uv#Az_ z3ri)`{OK9__la5|q6*Xbl6yz-gHo-6i9Zi+)v&dI#Z2MTp%?c*j-5=XRy~Wi;f~hT@H5ZcYCwd7H=@a_&?m;V!Mi zA0OHJ*MRAOtc{|*2B%TNJFM@9dvWQ`9PZT+axKh`^H~z{FFyiOC~Y)>L@g)!s{N}h zLf z%}B_6z<0el2AkA6y6NK!cm51Z`6aVTmy5EMcEH>Nbvimp8K1K3|L&w|Zy+GpaQQgl zJzu%%>^{RVr@sK+aW5{AU6uNB|FAT$(8m!px%j*)iSMX7<_WY9_OJIrXe|6v!jy6aP)$x%Zx&cHQQe*0+@GJIHMt#A!urc;bp#QJ8Mj zE1Kt=xtrT!Ia`}V$8!tu8r2-$pXeC*TI4;d!f#6fiiC_VuB@)(8=aE}Nv&+Q?vlv( z%%3(p&^xaL-^7%gHdo*ug$WwQSiWTpX6>VNE2D;9BZuA2{(2{P)dCRKMB6?_U-@eL z+MV%pzY}o1dI}Pw^LU9&ERdd1FI>NkotKiW3)&Z+O0cg-Q}D4kF5db&64q?6)) zy@!-R%9hU+=dShZ*=WCt`kVJjwrnPpS;Ot=#iPH-Q@4+0jje{4Nu=5(v*rBI=ix%w zdXpo@@cp#|4YZG+QW20k=#Vda;Kycgw+L&#l03{zs8b0$|JPTs>vu`Zl7}`B!tafe zLQ13YAAXC8C=>yFi{$KqlxW$ zy867orW}_fW{(tfIqhjkb>*om@laeS<=p2*)<<|a|6;oV= z)c%j60OD}Ar9~WUTdmJ$MN`?>{-J0Irlinr?hf1rYU$k1*9hXv#2Z@9^Zw&>sSDSc ze8E?)D8>RB6eqnrlIDr&y>k1vQC>)TW48vhN|`CqO-?LayNOJxGqeGtLVvr8s= zwyXjQ?@jLq-m#gKORV%0cRCl=Ag6`I=aKb8YG1c6GF73|XtfNnG@hn(bPC=neBQF@ zKrTzuj<}zp`_d;sm2N2y2J6tL-we>`PAZ!4kva{>T9$zr&(pV{iHC2{rg6IcBMx>- z$!+11pPela z;oCVuYEaqmFaZU6gJml|rWtjb3_|gNQ<|mFgHLM-W{Z=9zc`$s+Pts^4l~v8SwIJ@ zu1NNwNv;arhLQ*5otuy9pG1juBCqW|2oJcAyJuVFXe|6JFrHKF!N`oF|sjwIFGBkEq@1|SLB7oKw?d(5X zMp|H@ry}Gjz7Wyav~O88r)$O~Jiu!+uTSs-HYEECJt1BnsH?r>_7Z@`tS=>WUc;Gg zM=vc)O$cOa5kP6jG!6ShwKmjDGn#1mpj+3U)Q{I~frpns;d2?!LCA|%iS79|(zVw4 zh0IR&CyqfDb1F#?&7N4FRI5(;?Rs8R_xXnbP1oW=#ep)7T|Ren{un3+V3bh*L&!|J z18!VD%B>9jejG@vXuVXSP0hn+K@2SYWpA>;?BB!BDQaA)x?N~*fY7YVHI)awRpo8W zc7*pd-$!#PLT})Ta_E~rj2jHCUPye?o-xB5bx0c}D@+f{K!@)D%};MjIA2wZ<8*J^ zp3U?9FgCxnuQO%XGhF@}0I59MJYzme`P=c=8d8@|3=m3IlPptjsd9u&``J#PIQRS3q8Gg{j2ZBTF6@oY2KlSus|OVnqDYgZ|C>)O{|k%J zkyHxn%&X5;1VHT_Pe#1~Xy5yf0iR#cUvV`}HB-8j;6i&MfsbKye3Z&MwcQUo!Pp|M znq;NPbst4b37=Y1dHRChcskIYP|L@Dv!jEacfc}Q9P%J>+p%~5M7h)XhT8zWd@AMp zS|?e&%o)kKxo*Hx!%XNhL9|UZDYnS_Vuj7UzxvO0HlQ0-a#kX-`5WOcHzCIg`$J5< zZ8YF!Xya6pu08vsWmBQNZ&1ah{BwZB^u1LY&W;~^0#9zXUy=QUKNGI6>g;X4H292U zX#gRqYesDF>%GX>DWsahOYw_+C6w4ApH}lX5YLVvM|cIp@S9~~=xnrOUc&M)?=|GmwG+cx!0ctvwu`g7;C3*#lxoAtsCVVyt4*NPQa#5zvscCt?) zPw#(f^7LKq=yxq=J9*g;dct_~s*TX(CuH+xL+QnE4q{+d!DWW;5(btJ)W(|$0Q;)wTe+*uoyQOYG2lz>Es1tv&b`amm!S=Q{}s<%iP@S@R! zSj6gWNOxQ^wuFBs=Mowyvh(QeKlu_Iw&dotk^&0=2 zq0<*kiQAsbF|9s7(b#`IY9Il`)PJX6ms8GMvh@B~PvWqH?rlApbPImd~=ez|0k#r9cSUfxv2 zI(Wy5?|Av+3@HO3OKMXQ0Mu_kQccUb(-B0A6gj` zFAAy_QVWg~AWNI4w_ulSVu95C(v+m-1N)9j7g>dt(tD$}YI+)B3=mi>aT|Wktt)U;BB80?qOLH>Nlh3C z)hh^|^Rw`EjPO;br{I?Dpadd*`3iHD;TY-MGwCi2x5mwFW-No%INaZ@c-R$52f#iQ z+T`Hg2Ei31A;UPU;rjC3=|g_^LT*McWRkEdi(&7ia$;b9lSnh?vKWmuWDEBE$$aBZ z;>Rlw)&nLYd6mHAR-pN#t>wsl#4A9ctnU;^e!m-_Tz@N3A@6&%y}}hnJj~?+{DSjI zr|9^uGp*!DSW)8@|8=+y61)wDh*eQyt)8n~GC|DUzJQF%dgC;Fr)Yii_yB;}kiS+m zBk#P&TC~sF^r0o_4eX*5BD>t5YcDW{Mhdds`PO38#TYLHqG|MbpA++33ZRH?mwF`e zWBLUYb<%=Te6=^i=+G#0j7%ucQNO5gwg>XjSk0>K>a&(~S2-_#gCcpa^AB6qaXe;E zx^!gYZTP6ab7=OTTh`OJ>Gx~2f%MJA3wXzU7*0EMlukiB zD0j)9Vor8Y_imPWf1}h0>KOKtltQ8R)qRRX)Za9bxJi}G&C+2PH^02_ z?OVzqI2^ROfwSFL_aAp+vRAB*_`Tm_+-0H?!7;6nvx0vNl40apA0!fn3;GVFaAMxQ z^K_8aq#9h8Uwz#_4n+{ z6`X=VtHuIN_fE#naVsESRCISS?YEQOqeE3$oa|^v`b0q^tR7A$QdZ|}b@A7kP>8@6 zfHe@T#;*mu;ayp-o?*N?l0e zl~kw#@CKootC@D{i0M^m@54KYn40`vre*kcse~Sq6!e+Ce|tW@bT!;wCP#2S-*>S^ zZ_)@~Q6`F^3<`g|diy8%>2ie|@Z-MAztlOUr(%n%*p{=4`=e%D5}?!ppk%6!MmW=z z=kz2B?;ln>ry|M+e9FBwVx?Eb#67Yduo15hoMvw+*y44%Ys6FHbrlFN!dy>(zm_KS zyfK&u_Dxfrvj)7INNa7p2ku~17q?UHt;0NCzOBQRF&lP}g8n~me+TY>DNRRk_1^y1 z?Q@ByQb=7M7Wn^2_)pYYPuo<9X*}$a+Y`lhE5pRt(k8DbU)@ELrpR3coq+7Y*DOYH z28L!#@@luBVaU%pceNc^6um|zSYua}T<7=Gu#%TdpsCLm8+U`pS?N~avg0%N*f{Fp ze`P_U({5#+8$o)Xog|lTR9=g4$msk-VKK3%v*f5jd?ceKm}?7-^edqS7BAM}Xk_%3WqApZ#?%SWXH<-fMa zLw7DIx2z)6Nrd;PRSC@Af);0(q$BQpmi~R!B3W2cK_v;5A?2T`EorYK^#au!0{ z^7?r5+Clu~Paw;{H{ucMaS9T*n-~;wL#oX~aBMj8t6J?olVBDAw?-RoUm5%RF=RRO zw*Lb~-n+n;8x6)CjOt{<^6b4n8ojwjdmw$j$iw`=fWFNxTlMYV#d$fuie0iR;zvXE z?&@lq7)5XMyF%$vw9cs-+gqTZiElf|FS16?m*r%>DcVRAHNnD-9RA_&85CMI)EB0X%u8!n|s67h1xrki$#Zlu}$Ub`v$e`7$q@AUqI8J{Y- zB%WuAfMI5{2>ep6L7}KBTC!i55;_-C@dPnhpW##14F0P;J&U}+-eJEPg{#Y`U?1k?vbcbi>+>>XgN=wpVTVX2crdqzc$>H}-j=2k5G zSvZSD4<9l-H6$cmFejXmF4vwV{u}kT1HNm57jAy;gwvh1P0!_U{m(KBI(0}JT1i{^ z&sV-?MD}gyC;5gFXq9*TX_e`jcr9HcpyaRi%A&ILen3L!cD?SGOzxSy#YT{7EPP}E zHJ+xrov@s4urY`S=3O_xMcxt=AOHNCaLlV8fMuJ(DpKH@NgSWK-^L;)S{#{kz}98h zrY`8NU75EsGG8R`%#C$sG&h^Zi>~}KZk5ks+`GFE5EajkIyR1@kE=KqOOE_kus*7D zdg_eAq3vf4S%;i-xwq`KCUGA81ZqHLb?F;9bptDBglpm~=wkb(su++IcFE0#yT&mwX$oDH` zf0R^Nx_E{Oa_#-EW-2%@aiKD;O4h0JrLD-@i&@8UY^|ERBu?7-)c_zu!qM=Xex z=O#sdV>g;BJrC14x6{R?ZPEC*%=qTvRQcBSyiYpst1OmBrRz`Rsgy5!9);kTKi@zv zN{)-P*YS{!iAcll?i}lwYZ$=cNv2m{7s2~T;uZ46tM_lLwr%&0I`R9zKTNLWiWLEO delta 16821 zcmXV12{cq;|0YBTkyQ5V*-~UTA`~InHAW>#ma*^VDtmUaFOyw_?7K+@45H9@45FpzvX#;zxQ(j`DQ$M+@sru9!@-ZP){BnDITdu3J)GB zNXha%eDLr=ya8YQo2!-XAPPWAgN)2Ob2se?bo{~1z^tD;3j8J4cssA1K5qDai)m%L zLv=M!&iPGd=u}HhGIO(Q=8swK<%+{cp7V~=x>bU_=3rTSzVHI+fD!~6%Ym>(*Tyn6 zT{`BKw_w+7@{0F*sst*84vDG|mR*XT&xF?M_X#L;>Ht#p4f*ErwI;wmpa3n>G!9)q z!CJSPl#rOpF1%=WxAHHzpKqcCmttC{aXe{!J|4$ zAd%|dy7c#j7fBoZ}6imLy4Jpf5Cfld>^IZKgrSHuK&6!*utWzDf^Ol!( zaIc-OH@$CTG~&rZPXWk7?ZO;@oKT977@*?X93sK=ViOQ7T7z-8V@0ZmS=W@?K60%+ zI!UAGK7Q~vv*T&6$1XGD;AFhcEpaTPv!Kgv&~JuY4cGj~H)i^sC#w*B-W^kXd54;z ztF4?2+E27v#PwmfsZpk{V)Y7QXC=YJK3fUwqH)??lT|XD+IdWO)AY77GtF8cJ1ghv z?5Z|N%qGP+2URnCq&p)OzS!aGiSy?Cwg69eW|)>i>@_St##X;><1#8As0gcDT$pU0 z8Ui`C;xWG1CKrD1Vq@gA2hlO#9C);xQ0H_x*9jmV?$>nlQujx85xPRYd)kmr^>V#48M~5&4jfkWwMAlpLyo_vp?d+xgqM51!1o0f)F9m z8Nhx3!q=ZuuN{yN@!PRKnP=OqY_8PJ?E0JG!SnP_f=z59W-e{!aPX%W71>BWT1lC0 zQ=aaRwL}8lC^Tdf_laMza7+ElZ`xMVM+cA8&+jbB7sn3V^&+Q_9eD4(tN>ejVy^V$ zv2>i>*ClQjZ$U4PxzN;;Q>e|7shm=x%Nqdx1T*4~*?N5PIgWcnj!O{Gyt+j$Gvlu1x}HXi9a%$ zr=6ABU2x|;459ZfA*jbUq82~X#Lc>rEdXZo~;JhJW}3aCKP(AVO`a8yq}#o-&|>RTbQjPDjp)O z8IbR*QML?V>FdOe{MCT8jz_qxAfz3jfj3^b63>0UQ|wX5Jdt)A&I;9+Z%qZ4iF}X8+Pv1Z2MCoiR4%qzZd!qpDCErWrA8}2j z6@WiJ)WS6g?8M4{I>g8h0Iz&WNj2CaG6|J#b5I+c-dTge^cmP1(do0t8s5fj-YfJw zu&}hc;6}KFDWPB`5b^s^DBuRkYQt}(=+a9}T&OL9#cC*`V!acJ?woR4VV6IAi<-de zw>)UOcRMEKr_yN@du-pb!cK3{Yj>%kGIg6HNUs;YhlSF!jDVS206HV|ROOf{pfw-Q z6oTJN~Z*+sM(VfYNmbbCQ=?A_iB}$34xlPZmCr(F8#d9r;4t5DVLkPE07-3QhkI!EX6s?G+KU~Cj=6Gp*UF@MW zbZa%##hg6BJmMEA?MyBB)p0p2-C9lEnxTC&DYaL8L$2fZ(6E-)xS`86`4=I|75DjL zr@jDFzd_^PtXeK0`i?hBO|lgOqPkjLQ`=H7TR1!zruEpgiTp9YM_LfkfYe;bp zzejaZjq=P2nzGi^-S*7g(=W4A6QpVjYWZdr`CP+c2V`@ic)(m;t5Uy3!BgF3Ub$|b z+^3&o&dF#`;9a3-Usu*VZMh%X{3VH%`ZeUu6_MeGbGXMg!DOtO=%UiY`-QBys{}#( zZ;^_+1G7)0QSVZB((D3*6yt zhDHAwiDB>7J9aPyY%l!hpyZ0-8*zMn!Vv;KI z)GnI%t+u?mcLRQ5+ZU;!$L-A<>?dIqY$SU_)bY#x$DSFsi_{ouH1NL{{J2~CsI~ky zgP^U9XDSV!u`1_?llDXX|KTUKU4QN$ek_GK)qj5B9EnkOwOg;*xhgwSQe7<9wlFV!vUMX|UxbD3M<^Z34OxX3;X5kH zyQ->jqZ=sR`t9Lj)FMvd2QG9q3!*wb_D=)hCOk9djg#4|S%H7cH7`NUxVhcW&ydFv$1am6V=#mGA$Aj*& z<)zE#mL=|I*Eu;!OS`vH`7^s_UiF+dOF6d8+rmeVRs42{ShwD^AwRC$JL7H~R4qaif4Md)pr`DB2^87cGs3ORXC+ z4qKR(`gPQJdF{4+E28BKO03E*TI~)LK2WuiyTo5IONGprUn>5auHxCfqTT?~RaHU5 z0lc*(=U=W@0isQrd?K{YTkP{c&XRVHJQRVVfE0fn(oPk~GD|TAvNzfsw;H(I_?m+bmZB&ghz`%?*^jRj z$HHXgs|2u`%LPNi3_rq{7e$tK2h6$ZYe=p;ull2$*Q1Z5?*%90^IVAm|A3jfopC5Tv zKHd-u85q76h`YGx;?=5qod8WxR*}NF4BBLwWgKU_eyr z;zdI#OhZ^TVCcewO_Gm)OmTUN?!2B>a&J9!RnD{CEiFPPewbQAYc}IsCf^?^fd|h# zUhH)D{U!Q5rf@>fSYr`;5&{ZX`J!{7SQUCdb^p&_o*6|zIJP%i4Bhh=QV!q6hfkfW zyvPR8Qup1JTK|EGWG)8Sv}!?2t%;K@<2;~>n$fN}O)2$4r5@g+Uzm)2`+&Ql_pEA< z0!)?*bXN{;Rw|RIgLuE5tbl%OKEu=a5$IDvh^cwf;>9<S;|dlER3; zc(EKLKV!nq$A3TBd*QPD8M-x>Z*O70d)O26(QD^>B3L$n)HJ!}fZ2Cvc$I~YvESC3 zLGHFa@Cf6c2`HA>9y1+?O0-f+E{&)L-Gvoe4{j}>)aMnv&vAo+gyVS6PV@1+L%YCR zu?QoDSp|SRaqfm5{}r+~1}93J^&!uOBVdwE{dLTJgJ&p4R!#)#>AWwqydLh7}=rAlXn`iMH+GHUYqr4yW()Mo&vEvS_VGPKK-);#jjlqvaPfiksR-d?E zf3CWIW^K8&>_Vd~SBS%gno$EW(~^d-fMo8VeNdDh-rn`%Ha|FPB z*la$p<>OOQ*!|iga+ zO}yXcp?)F+E*hGIyI$P+J?m1L`p<7uy=2wD5z63dIng8|cv~HQAMOQgENDGAUxCp( z(5p8RrYO0%toXtfubiIls7p5{{ggHLd1D_QL6J0k9?yEi)0*;YAGO_hyh;htsZN9u4lU2i+dDPR=e zP@hxWRWB-Fr27BlG!^;ZFu^GtrF+HY$;vbkD16Nhpsa^9osLOWQfS>wPkAR&XRu<# zT11p5q($1Pv5)tc2i53S3zP;Js9=;qq>KB?H4fCqppquT(~b+%n(^7I4JV*e8Z-I;!@Ue&TIm z-jFYSO^0CA$5suekDa(|+s=^B3dq8;DP!p@Z@oyjI8;2lXoH)G{#(-wxV3`l#kWJ1 z!qio{4rE;F+VY(ZR1*9B2IN>nh;n+h0|UU8_nZhm_i-_SarEg|V>4^9x@U9xTyC3q z&VRfQn2BdpJPL~OugS!=?i~vd_Hxerud^OL&N?3Ocz6&L358;Reu8H`hYef6oOt!j z`kyfMfrkIPv5^F!KEVK*5flx6?D-SCu8>X;N@<8JkHb?Tfn(uSTyHiAd|&i&VxZ z8hu>R+HK16T;QI;%s5WX=&Bjk3zXB{loB%^vcG7>w%P*t))XU(Q8TjfKxWE~nN$16 zwlkU5kJwP~2pjR!JInX(WM?sh@?WxY^Dxf#y1YMNq*O{_?>y#VP#~1vIrkw)+Hm*^ z#lf6Ci`SoUKkIs{InuF?>}9Y}WdMmTHMUS2RI~OjQ+ItX%h7U}Zr{0c4rWxLzH(RX z&b&LL)7?!#L%k1Dc~p}B4msknek!u{!G$95MO)vq)z6fAIzu!?jMM^(}`$7=zG>uMqnln&gAy&^pZ3Y#zNNP-Q|4>_myNw?6?SD^|J*#)iAHOUUKziYpll0JSDWxwz(s4G0KD(%b@ z18|EDw@IDlKLez!b~IXl<`z{jj!Y1Q*=LABVooN%1^C(pFah1|#*e&JuD0Ss24JeX z8NoslU}fLRAEIzPmvVyteo6Kdj7tL-hi^PH&lERvoO9mPJPrCa` zFn{#wR0jTnC$@A$3lHb*vqqAEgiOdR*QkGo4*$@=L04kwh+AkbadpJDU#|k%Abir^ z=~F{1cz0F!th|{qQyP}wqsxkAOcdq11F#BkQu%7(v~mNxVGY)WE?)~e%Q>hkukHC( z{A%r6d=p#Mlgi0Fmi~)!Jm%ipzUO+|nB&94N~DVdrNqFm;yfOakWH043qntcV9M{e zYej_{TnwmZ7YEJ+pUpCg?9q)j=fin_85BKTF?r~m_!U<5_`BEqA9U~Op{YBEw-4Z* zEv)p)>7{FEYOKif?#j$|LnH6}s20J_r`-CSJz7d|PHM%2pg@DnQk22SOvSfRpU4Yx z^(OBa<@eS-eL>s?2Yfk1WyTu}M}w%~FDVbNwXk%aW$9j3UPD)te*NKF*Pq@fQUxyS zZa8Hh`JYQGhRdH6N7qc=uq7WmWOoA=`zNe&C``$2jT{>eNc)9G>ddN4mL4ty&kg@> zH9GlLWk%pgEQB)kUw`qSXM|NTqbx7%#}A;wr3A&St%gt*=WqBNbs0(}0RaqinJag@ zBq%dhaP&>&UuUuHTuA1ZC*UosewlZMCKjY$U`quE(X2$wXu#Az zPtqx7=Lma7+>$B>b)WV$EM_&7TZdh>p_f@}TYkf~KF;T;qTHEtUHhf$pv7er?PtRlg(BeR#UJhoKNziK5J`%nCr<=!VO3VR)P3> z!~Yc6*h%W~Z>>1mcS7$pH6+o4EeLaRdHZ)t%+|-qZRA&-EjY08qJKVm9VlETMCEZc zaMXjNrjvWcnsfA%Eys^mz3_^m(dND6s%+UoFZ7P3kU>DNmQ6BP8L!EFXq1l73N?GmpGyIEL^DuFf5USjH_71xFL zJH0u9x)S4#2e-Pqs%~}Y{3&{bf)4qN4@T%Kf{aRz=~!*R^)ajK{wQ4ea_3TWi0J)> zhR%lbW~8OF!sR(Ig#Ni)`yqC5aMxhtiPwkyuTX&?Nk=PIcN-VaEah~HMxU`W_CD+J zrB00$BgS}i&a9R!bg$lJzFxWMGdp`jxtE`cB2V_p{Efhj)v%`6{6m7t!9omo{RCR3 zflIc^*?JTY2R+3ZTQx@hE4-|DGrG8JeV9ML+(wZ=H9b#caU}JYiay*mglNMw`nRSY z-NPBdCd)&T-FWHugP!@2x-$izKX1a7&L5S9JX;V_q{)$eC-?D`!}g3*h+|$tff<8- zC*37?GcTK15kE5wJ!}=YEUfERq=g(!(4y8!)OP{yo50wmnp*gC_Ijwa1m zQPN69v`D*n^sV0D0)k0C9~>?L2XsU9XRlE{*!XbhZ*(wrHXbI-!fJDjb!UToJrB{;{m}%pe!Ni6P=+K%v@va4;56x+~l&0HLnSjyJ26K9X53N*|`qv`dTSNk7elt;G-nG%=mJ! zI3A^LbSCQPonFl!-F^#cN5i6!H9J`;MW-I#FbV^`Iv9%?J0M?yj#pF)(?pnGveU-a zAjYo&PtSH_uPFGfeK^uL?Qp}t$rXq*fQm)|F&DmW%ew}8&w1|8bqX>dA6`1TP+Z?HdkD2-WD_$}_1m;X z2_R#y#LjvhwHoXyo?TpnMxdf$?7LeQpNxCCEz_Im1DG+qO^$ggq)y2xDK?dg5 zA=|Z1z*^$yU8KVh;}U+rvMe;677en+SMafHIm_4SLiHJ}tVp(2Yqt%tV(rj+WW}ymE z)Y#do6NOX;Dqe?|p1LBE5o|kp6BXEMD-z{@p_=ykzeU)EpW>Y1c4}96rE`>G!i-%~ zp1@xk+_}$a4Gx~B80)=TRQqi@Y4rsw@P)gJ4hZ+!eZ19GBj@3@Wu|Q2*LRVX>Gp|x zu82dcx97<(Dzf58$gS3UzBlM!izLzIqMlb&0ggqk8JailaXf3wS+2hc>F19aW>REW zd?2ZIc#t+6laQ;`Z}kh=%)@PE;w@6owD0bt0yE*b-F4$M-i2fKo7CTp)C97ykLQD< zn=`3n)#O*^5v-Z}_#c5==i2gW^T*2jmXyF8q1Eni%q2vm{K0vP-J(m$$PbEY+_n*8 z6tEi-iV2FzTd#V&CDe;kF;qmmh+oH=RT}n7t#BmjIk5CQ{tF+Fa0DNa=2H0y?DAFg zirO`c&`0Ky);GVd_xraJ;pd^Dh|+^sm)?ZQGb#4XINfcZiqYezqkQYhY!M(Xp_Ppj zkK>8a#kUBp)!>}UMdMgxKuC(~gK5SXK;50?EF^_fi|okK!{<2Xh)bjyFEdxi2%f6r z4PN0=kj$36fek3ye(wXW)lVt&5T#U25~|fkA{cZZT8fTk2g}QZ|c=n^NS4~^&0Z@6!jQ{TQm6I zJW&z^B|)53-Gq>VB-Vp%Y>|T)u(w`RYcxNQIb-sJiH_e#n@;NsqEJ9vCh(}Yt zJ8A8&0z=+(gGBh)c#a5i3c)7kFY%WM`L(IEBGzY&!8I?^n3gYo?6+|oCDca&$8rlW zC(9zE?dXuMGlLm|ACJO~amah#3ch`0>x1V5;Krjq<#oTfqhi>Hj(C8}lG90Rw3_0> z#B%zdR~xPD)AmkQENW98XP;^}(D4BXt|6`k%a%_CNZa^-3X? zQcclloZ&S;Lm;KGwv2FnRrGbi96fc1R;fg}vi~GyVZ5E(ZF3*$uQ}?L^6yKA-tBYb zT00){iy>wMt9Qu)*>dlE_j7n!ryIb1gJH7EN4Mt7N73B%yq#pe*RDIXl8($5Q{fW} zNcs!m&C9a%@6|v5iGD(!-gNn^bpLct!{7|@>=^Dz#;A92g=fv&4*j>Q>5n7MGBYgNQ*4U(j=ClLWRIo^lH8Z^LJ}yQ$d`Bf@y*C>s8| z!g_3RMy7^e*YEYj)$*)RD?3@BvQJsBy)ndQL~})YG)xxmjD2%~tHWyz4rZpm3+b`wrBow?O?J-te#4 z#OoA@oPW(ZBJvspIo~1PClL$hlJ>se2;^b8bn|c&`Jb%vU#nWtLH?n^A-z5vvr%q; zZ>5gX^z;tJ-t#}Sa+ZMU#Rc--bZ$qffS|CA>n;!sjJ;f`GYNR|%d)K51?dUi{b@1# zBUVR-CJGEG2Wh>T$u7-WZJ?mQzWh%=%5AzAiV^y02U1;EI&NF#{CWm7*y&;U`i0X& z-6ArseYKf(jcq0|qaGaW+^tyi>Is8|`{;+=*G9e$%GOHf5sCbm;JFkdvPhSbnal zx8Zw*tnX!WGV$`X>!N~t$^5Iw#J zPEOv~nX>-AL3LBOJA#b&>oAWg_WcAmalzix=PI!r?Ims#hKgM+EjZjus#eh(TWrC+ zE?kHD0-Viu`vd*nz)Sn)ea6In1UB||t$ER78T{Xg(6u#stM78{rknf%QMXQ( z!Oid7kyzn>r^bGeF2+DGOY2LK3EHX&MLFk*mxYFlo)IRJ1_ttXP|`&9p0uOhuEsLY z%Wd6ct@UI5b(E2wUn6aM3<1Wldea8mxDS{v=GUn66kMVQZhuM3<~}R!Fr7Rrv-3?N zv`x&@meg7hEkHrkWQ42#%gX*+)ac>te{EV;TFB|C#!=rPSX?v4H?}yPGn@Z>&S_9U z|M#nRd~;hPDh8kxr;O|G`MIRKYs5;J--ietVlddAjhq;!PbF(PUB<6xdw}~guaLD1 z20z;Ttz?e;25g^Sve>k?zAy`|b%>-r5$@H$-5no2Orfqd$ER5lNcU&6 zgcB{V;PG$F6fKpruS9o+e5{x}O(IQ_@&YhV( zch|->?O47lF=TCaM#B>Z?BC*|x`_{>*Ek~|+f?Za$5V!F^Dm>;{efv&+w6MbykB(a zDoI+JiMqA6iev(D886?mhm+7p_m0S0Y|?9qMe*I=Kw4r&$~tBk4T3LU%vSCw?=f!M zZ^xg3SWmBWG?d&srAs*~7w^|8XSJPb%VZ4jz}|)T)re6k!tvEkypYC|wIB&%db+DRB$yK+>$mI52 zIPoe&e+9?UMEKJbrtT;W%vlkimp~vfZT3nV)#X}#hkj+|Ghm3vOt)jd=!LwaDGasa zkDA9`Y~a2Jz~#C78@%tZp(aV+Tn>TNmAP`hKI>BNF#Akkc);0<9qP88VE>v_uY-T= z62ws+HGaoM0ov?lRCUQ!`gEeo6had_ffeQG{xxX3s?T+2_4}S|1qj>%WpvJ0dDZEb zr7y1_u3nV1*|HQEuq4`emh!&B<+OgG8dBjkt2|CPO&r2k!>lf$jtCHm{J%g>ef^IK z-weE-c=$jTcuG5LN92UB4#?LXlLcP!d$hX{t9{9#@z7;~=@l8nO#K~V^b}I#^oWiC z`GTS#e_}YC;7k?gkTT$`#J?;Y*rcSP;vTnCa)71G9#CGvGmZZsI$i^W55V1gq$^j^ z`*Xj~-qL)#;Fa85J~CY!A%D4=5{Lw*k4bFPtZj8&W3Cm{LwWDichMb#Z=?!3Cd49(kV-B-f`PITb2QCU?Iaz1^V~j&>~DWveipTwBG}0+snFQ&Z+enu zUUIIigw3&R#t57$x@ob+c5>5YkRpDLt>HYP(cvpl>g*ViHE^Mv|73wTkB%`}VLgjW zK$)B0H-<{NC!WjmgIMLmz>Zhwd^3^ko6+wY0j1LV2-x<1^W&jd4{Vr3kwyPm79$n$ z-;$wkrDquFTz75ny4maz#3dr|9b$Oj^7i%HIKI3dqS4X{T61-uD4)VAH#)3TeZrFf zCs#o_&??`2wWR9xk-A0C2)xBDHk}l|In_#i(Yg2FgHLIkE>w~_0q4PTCa|-CQIG)U zWw4H*0F?tXk>N07qz817d6rICXW$5O$R=O$AbmH8obs#{m_Vi$$ zG0TUX*cSau1<@RfTIlJ)9+sB#-HrOTZ`|o)Z9-D2m|qJ>n}+Bfs=N$~gd4A-JP10g z@9|X&0pFduQ_UT@St}4)yIl@5)qwDO>RK>Y*KMUQ0R?FUHqOq@DWiAs($>O)=c#a z?smyazqVfQ@x^pyfJ1!*Q_0Rx8Y(ztds-XfZ)37j#dkSlmgfB9#I;Ya9#tyqi986Q zM-40?#TKW;Lo3k>)DJm+_!Fu;nKEDQ2e~PF zgJ}+KE{pN-b|Mw-0185M3PQp>iC_FE!6u~Aysk@Y2U=RZ*|POgY+05L9oX)oO~)rI zvYf^1aKAMXZm9)o&e9W1Da2vIU6aw<==&3i`-W|?% zQor>_WDVCG1dR@wWd?1a;Z<$+a?Y%}>D=J{<^cmp`Uh~^McJzFMyU>`hdNk5Hr1X6y#xzRhNHb{;_*R zbfS)s({r5D2&$bRwhu0Q0@LpB6WrHorNcZsoNZOu)|(fCOxmqt{k4rZaBe6OE&LNk70L7&RM=mu-v1EydecdRmPd(I?~)v}yg+D{vYMZn0Z z-P1u8WFNRudhcD9X=VlIng3o(&5Xi*PC(NAgYWc?ISu0}QJjJA+2X(>xKe#vP$ejU z;CbpX(qlvtxufd=d%syDmSs|MtVwqEq=tKK7BObFIyTIm;CpIS9y73k#o-FV2QVbFc5^6a2F-Q`eMvO!P4@dJ_qKm$ z5M@j-PWJAXLiWumyn2nTx=g-Tp8AK?ATK^vXEwhrA}ez25uI0jkRR`s$QsSeXwr*a^`3mgefXD-GB3UDfnDGnL zVL!^kNStt4##sKAI7DsRZ3zZzhEbyxQe!9ITATV@45#Wx?kv=aoixpIGIRnq+72)@ zO2xyniJvVen}djKfz7Ql5u(9=yIUH31c@=?0*aPCPLWmpvM!;luVq)=YBLSv8H z)sp8MoUQ!cazgUbSK!Xs&=VN<3}+|yg_fCRlN_G zjJ7o-=DkLiWpI|Bg-PgnfZs!skE1j7BMIYrp^)mMo6;@-Te!|F`1q(ps$6*az4Cy+ zW+LH`IijKF&bi#XUaAUc*hs?{fJCh}Y#$|#T---_2x1D=CMad!$AdSsWh2q%kNluh`0^*)*4 z8-?#3!t}tZBehbaJ%-3gRo?Q#i2{Oid4MEMp89u*>}4|Wu1^!LC+7)Su*e?9j~UY2 za<;8Hee`$>;E)?)^k<;h{xYfXZkfz|rTeGFfU;NG#vkYB88W1Wt~N0+gtNXigeXw{Fz{1bRQ;fwp@a{<08_?j4Q zj=NuV{JDBQk)}xhPf!>BV-~-&#*)e3rJu)ntHFZ9LJId(;Qent`UO6xr-Z&v1}o@9 z7z*s)CIk^g33sOLs$>eyIlQNnHTZR}3!gmMhvvuCl0hU?x8EGWxciR$@UwC@a7C~mf3#t*_;}T% z^HG|A#Y#nJeo_3U`D=51T*=w%uVC?>iIHwi`IK50ziJmqQkKTN-ZiXRm0_>c1I>@1 zZt3@_FaLz^CH^xo7NP3y2{$n3XdLiI8f;E@EC7;a0pXV()Se@D{40V1;X}2{mikCb zW&wq#h7s7J1TJ32%!A*R6ydw;d*e@TGxEue%kIE@#17TY?y01n&&xF7PT5^%4sXJy z$;Uw1@#kvK_RKag8e1xA`!!M!^Dgh&NSC5A3`TFj!3L*OSNh4a-aBXP2LLU#Y;8)7 z_RF-_22VaLH+j%@A6(Hn{2Ze-m|Mrr5bN)KO03PzevZ@H%;|WOzm~q2CZk+G0ctvK zFlX=4#jB?nDkn$l%Sc*%vt7fzkA}C3{VFQSMx}%G6FX?ZlM*uAAUX6WW}&&6CH;p3 z?0=9_)pySMKlD;T-)OC*11HXfB=13qyr&3l&bHVi7rwQCx72maupASAghKyRyLl>l zUK&`ytu{4wd~i_;&%EEGK5@Es))QB<)-2GcGK#yq@kUjzJr?(g_~!QJq@z&25K|X_ zK6-T|f{y0man6PRrjL|x=R`gZdXsrCK=m2h=O!gS@`lc1_v2U=;41mcLQCv<<_=A+ zGDf|Mp@nk?@f@~w>AM&Tr7!wAFT0)PP!~= ztbUkqn>jbGu9pA&=z3VJ!tsn`X2>vly6)Eg_bU>JuQ2KR{Xa@k((^yAV4o%!jz~QS z&4V~-1e?qKq0b`&DD<@$fY*kXKC9_04TU%2X+9k8RK9D&GwY2&U((#%llSAkW^;$z z81_-epORCZOrn{y02;f+Uom0KkP#mrEdz$i^ zK|RB|calFNKaO%wPAl1llgn(fZCn>|nXJR*vD<62o|2udR&4Ivni(#GU!mKR`CBo)y1q0l-DXVPkw$%@z{{Q5XQb#~&-BySiK)`ag%e8gH^1 zzv5}&E&qFpsndTdAT7AbV@6<<>Eoc^WF?n^3;A0~M<0*FPONn?z+`GeeYuEC~K6e~w zWpnD0Dc>Bn%+1w$zVwe@JB%ORg*9k;OD`ri09RD^?hhnN2E!WW_sj}s8olk)#Y;Ki z8~zF~#IL(1>Gc{!lQOh;cYc_CMF`93{&Vl66}kK;1A?Pjk;8~qgSEZ^#>}4BmOFjP zr!%tI>Zh?z%8Rm!o_mazncE(Txa!P~dO3`uZZ_%C>yl`?Y^k3Nvhk<##A>l%y>v5y zVlRyQuwmT6w^a1F%jMmJ^UpUMAnx^+*Pp$;{AEw5>%rH#fcghkKciKpxJ)dkEEkxp1W1#VL~nLAjOUl4-s>~l z*vCiiJG*h`SYaBeL$v2+C)JTA6aE75UXgRKtB{4$(@ZYvtM`p*v*}lcIqb?c;cTBN zG2)DVFql#)rOa?_F;G#U16tY?qCAYVxu<1M!lP@t|23&@T|R@%n_aK<2#;_7f$@tV z%81)rH=LNskjku44>KLwJm#kc_IJ{{KOGml6Apb!kfF*0YR;Yem>l+;HDzGIr5Bqq zm{IpiMS{SuVrC6f5ReG2Bf&~a5J54Cr%n2aM;oz+zSXJS1_PfBZBiD!Hh)}FOj-OMa$sZ{s0%lyYLWPj#3wY)%E_>)$B8q(a}Q!_Hy>RzM>$mo2!=dzx1{xQzuygXY0aW?FW0+FZ%wf+ zL0lEmILXw$gJ`9>Z(P5KQymCF1&ZgXr~E}Yfvw4>Uiz;E*>oG9WS4U2MX&6B;r($H z@=ZumPN^DY!S#Fcv&+ALe5FT%^1PLbn6$Y?dl38Y5y;Yi#37t7#vA+;w&$jVw?p`j zOVO&)e+$Nx_W;jtrU-Fl&oC4+P(Ab%qigk3jtjCFr~3(*uL%L-3&2JI ze|O~3o=8#TOyT0Fv>|9zxyrkQMo**1$gK$RAYXttQ0skH>?7U^Ss}} z(0W%Lfq15C47~TgOj#}7Tjt#y5T--1%`stO5}5Q2!$*YY z$m8vkeYP+wj|K5k^aL&8)62q3Xd{{w_m`-w<>e)`{DaOr-N*;bXvT5n-HoL`fx?fS zyKj9UJwnVMq4M1)DG3@Elj=4u{%HYy?jJRSp%cgOsZ_k#hL15EDH(}!x15&VFHMIt z{xmtl0WUq3XV`fV&$1R*Pb%F$`j}eFQ8}raE|}&B&BccG!pS&P*MwlbYg7lC4(3A9 zS}`nvSA3FW7&~dSF{__(N4mH4zf(JG%2l1(hAT0i=f6ZOt{)5*y8W;Zl+5%Rz1V~O zoU2+g&U-k)93D8KW^iLve2m>7v8|HJwr%s?s6DKZK z1mFPM9-qoK;kAE?5lUt+#O5?zQC6Ind*J7se=69UmGAR(6_>gC_H6HRTHv?4Fi6v; z2S(jIor_~L(?f%nYx|I=6DRb7+#SW_NI+G!C+!2fd)9lp%HTCW|pw<%<`0mWZ? zH}aW&{}9bpQM=aq<8wXhc@_tZy1|>)B%*Z9Iq4k#Q0`IKqt6Q;>-v|XIEKlhJ6wrd z+ahxXc5~Wbix21<-*)vUs~h>9zW-KEI=Sf@e!Ih+!{V>)t-{o|o6FE%^Mb6Y=tO`R zdvR0zQALS|WukNK^f1@mN^!^%o|d3~O>!+3tv=BH&+89&P*WivMQHU>JUfpKxdAjN z$;154j@%UvJA_s5U<&8ZmBDOP?5bb>3_i6vDb`j}yA0;^Yz^kTN)uZl+leF-o+??{ zVEMt#>*dEeuwC1nJfXY-lua!}=>t9c;>~Z1T)>n1VHFFG^}ec5)GFWDpKE+Kw@{`B zH*s#gUy@g~Ot}=!3GW8ngp@Q~j1NY_=Mx+@8{Dy5y6imxbLwoSb(S_9WvmB7+0YQlFh#ep4xR&z(w*NsX^h;5H+L*^B2My^nJ>Vnv{`Zr8 z%Bw3G`M4fHInP$3-Rc9#m+a_RY>h4_Y$EE=!x{z%_M!taoMRfGgr+Hclgs}F;s720 zUu*OILBv4W9~RB`PS1Q=CjRbA`iVN?X~G0gkr(o;_%b58`|QZ4Hz)|PPqT68_yrzxY&AD8O44)Xm986|&NmM)%Qf*kC+y76b^6M`r!wGZz}cp>vM(iIBM z-FI=G^?&47MRwqGNPGurts@ph%6pR{zp)$5mEMQx+}r8m(z0m$TV{OoaH{jx_P$TL z@2f0!j!M^`IHyv+?0FP|WBz;ty(qaZ(q6|yIws;6c6aw!$6Uhz4o@<@`U)4`M-s1) VFJ8TWW3_F&chrgB{{^2ptB$AE*8>0m From cbb639d2851856a74551176cf25990c94ec16c1f Mon Sep 17 00:00:00 2001 From: Jonas Ries Date: Sat, 5 Dec 2020 16:30:52 +0100 Subject: [PATCH 2/4] implement proper weighting --- other/BP_Gauss_analysis.m | 68 +++++++------------ plugins/+Analyze/+sr3D/Biplane_getPSFxy.m | 61 +++++++++++++---- plugins/+Process/+Register/CombineChannels.m | 8 +-- .../+Fitters/MLE_global_spline.m | 4 +- 4 files changed, 79 insertions(+), 62 deletions(-) diff --git a/other/BP_Gauss_analysis.m b/other/BP_Gauss_analysis.m index 08f6db30..74837cc0 100644 --- a/other/BP_Gauss_analysis.m +++ b/other/BP_Gauss_analysis.m @@ -1,21 +1,26 @@ -pixs=100; +% pixs=100; dz=10; +zoff=750; %% BEADS: calibrate sx -figure(88); +pixs=g.locData.files.file(1).info.cam_pixelsize_um(1)*1000; +loc=g.locData.getloc({'PSFxnm','PSFynm','frame'},'layer',1,'position','roi'); +indg=loc.PSFxnm>0.1&loc.PSFynm>0.1; +sx=loc.PSFxnm(indg)/pixs; +sy=loc.PSFynm(indg)/pixs; -sx=g.locData.loc.PSFxnm/pixs; -sy=g.locData.loc.PSFynm/pixs; -frame=g.locData.loc.frame; -z=frame*dz; -plot(frame,sx,'.',frame,sy,'.') + +frame=loc.frame(indg); +z=frame*dz-zoff; +figure(88); +plot(frame,sx,'+',frame,sy,'.') ds=sx.^2-sy.^2; figure(89) plot(frame,ds,'.') -range=[45,95]; +range=[44,100]; inr=frame>range(1)&framemaxrange(2); +outofrange=dsmaxrangeds(2) | sx==0 | sy==0; g.locData.loc.znm(outofrange)=outz; +g.locData.regroup; -%% data: average x,y -trafo=g.locData.files.file.transformation; -inref=trafo.getRef(g.locData.loc.xnm,g.locData.loc.ynm); -cr=[g.locData.loc.xnm(inref),g.locData.loc.ynm(inref)]; -ct=[g.locData.loc.xnm(~inref),g.locData.loc.ynm(~inref)]; -ctt=trafo.transformToReference(2,ct/pixs)*pixs; -lr.x=g.locData.loc.xnm;lr.y=g.locData.loc.ynm;lr.frame=g.locData.loc.frame; -ft=g.locData.loc.frame(~inref); -locpt=g.locData.loc.locprecnm(~inref); -lt.x=ctt(:,1);lt.y=ctt(:,2);lt.frame=ft; -[iA,iB]=matchlocsall(lr,lt,0,0,500); - -normw=1./g.locData.loc.locprecnm(iA).^2+1./locpt(iB).^2; -xav=(lr.x(iA)./g.locData.loc.locprecnm(iA).^2+lt.x(iB)./locpt(iB).^2)./normw; -yav=(lr.y(iA)./g.locData.loc.locprecnm(iA).^2+lt.y(iB)./locpt(iB).^2)./normw; - -figure(91) -plot(lr.y(iA),lt.y(iB),'.') - - -% plot(g.locData.loc.xnm(inref),g.locData.loc.ynm(inref),'.',g.locData.loc.xnm(~inref),g.locData.loc.ynm(~inref),'.') -figure(92) -plot(g.locData.loc.xnm(inref),g.locData.loc.ynm(inref),'.',cn(:,1),cn(:,2),'.',xav,yav,'.') - - -ind=true(size(g.locData.loc.xnm)); -ind(iA)=false; -g.locData.removelocs(ind) -g.locData.loc.xnm=xav; -g.locData.loc.ynm=yav; -g.locData.regroup; +%% data: average x,y +%use combine channels plugin \ No newline at end of file diff --git a/plugins/+Analyze/+sr3D/Biplane_getPSFxy.m b/plugins/+Analyze/+sr3D/Biplane_getPSFxy.m index 574a4147..9f3abb23 100644 --- a/plugins/+Analyze/+sr3D/Biplane_getPSFxy.m +++ b/plugins/+Analyze/+sr3D/Biplane_getPSFxy.m @@ -17,9 +17,20 @@ out.error='selected transformation file does not have a valid transformation'; return end - file=obj.locData.files.file(p.dataselect.Value); + p.datapart.selection='all'; - loc=get2ClocPSF(obj.locData.loc,tload.transformation,file,p); + if p.allfiles + loc=obj.locData.loc; + keepold=true; + for k=1:length(obj.locData.files.file) + k + loc=get2ClocPSF(loc,tload.transformation,obj.locData.files.file(k),p,keepold); + end + else + file=obj.locData.files.file(p.dataselect.Value); + keepold=false; + loc=get2ClocPSF(obj.locData.loc,tload.transformation,file,p,keepold); + end obj.locData.loc=copyfields(obj.locData.loc,loc); obj.locData.regroup; obj.setPar('locFields',fieldnames(obj.locData.loc)) @@ -53,21 +64,43 @@ function loadbutton(obj,a,b) end -function loco=get2ClocPSF(loc,transform,file,p) +function loco=get2ClocPSF(loc,transform,file,p,keepold) - p.datapart.selection='all (T->R)'; -loct=apply_transform_locs(loc,transform,file,p); +p.datapart.selection='all (T->R)'; +% p.datapart.selection='target'; +infi=file.number==loc.filenumber; + +inref=transform.getRef(loc.xnm,loc.ynm); +locref=copystructReduce(loc,infi&inref); +locth=copystructReduce(loc,infi&~inref); +loct=apply_transform_locs(locth,transform,file,p); % only valid parts? -[iA,iB,uiA,uiB]=matchlocsall(renamefields(loc),renamefields(loct),0,0,1000); +[iA,iB,uiA,uiB]=matchlocsall(renamefields(locref),renamefields(loct),0,0,1000); length(iA)/(length(uiA)+length(uiB)) -loco.PSFxnm_old=loc.PSFxnm; -loco.PSFxnm=zeros(size(loc.xnm),'single'); -loco.PSFynm=zeros(size(loc.xnm),'single'); -loco.PSFxnm(iA)=loc.PSFxnm(iA); -loco.PSFynm(iA)=loc.PSFxnm(iB); -loco.PSFxnm(iB)=loc.PSFxnm(iA); -loco.PSFynm(iB)=loc.PSFxnm(iB); +inffr=find(infi&inref); +infft=find(infi&~inref); + + +if ~keepold + loco.PSFxnm_old=loc.PSFxnm; + loco.PSFxnm=zeros(size(loc.xnm),'single'); + loco.PSFynm=zeros(size(loc.xnm),'single'); +else + loco=loc; + if ~isfield(loco,'PSFxnm') + loco.PSFxnm=zeros(size(loc.xnm),'single'); + end + if ~isfield(loco,'PSFynm') + loco.PSFynm=zeros(size(loc.xnm),'single'); + end +end + + +loco.PSFxnm(inffr(iA))=loc.PSFxnm(inffr(iA)); +loco.PSFynm(inffr(iA))=loc.PSFxnm(infft(iB)); +loco.PSFxnm(infft(iB))=loc.PSFxnm(inffr(iA)); +loco.PSFynm(infft(iB))=loc.PSFxnm(infft(iB)); end function loco=renamefields(loci) @@ -84,6 +117,8 @@ function loadbutton(obj,a,b) pard.dataselect.object=struct('Style','popupmenu','String','File'); pard.dataselect.position=[2,1]; +pard.allfiles.object=struct('Style','checkbox','String','all files'); +pard.allfiles.position=[2,2]; % pard.datapart.object=struct('Style','popupmenu','String','all|left|right|top|bottom'); % pard.datapart.position=[3,1]; diff --git a/plugins/+Process/+Register/CombineChannels.m b/plugins/+Process/+Register/CombineChannels.m index ba5293e3..e607b7a3 100644 --- a/plugins/+Process/+Register/CombineChannels.m +++ b/plugins/+Process/+Register/CombineChannels.m @@ -45,13 +45,13 @@ unmatchedB=copystructReduce(loctt,uiB); fn=fieldnames(loc); if isfield(matchedA,'locprecnm') - wA=1./matchedA.locprecnm; - wB=1./matchedB.locprecnm; + wA=1./matchedA.locprecnm.^2; + wB=1./matchedB.locprecnm.^2; else photA=matchedA.phot; photB=matchedB.phot; - wA=sqrt(photA); - wB=sqrt(photB); + wA=(photA); + wB=(photB); end for k=1:length(fn) switch fn{k} diff --git a/plugins/+WorkflowModules/+Fitters/MLE_global_spline.m b/plugins/+WorkflowModules/+Fitters/MLE_global_spline.m index 2349761a..5d389676 100644 --- a/plugins/+WorkflowModules/+Fitters/MLE_global_spline.m +++ b/plugins/+WorkflowModules/+Fitters/MLE_global_spline.m @@ -253,8 +253,8 @@ function loadscmos_callback(a,b,obj) switch fitpar.mainchannel case 1 locs.([names{k} 'err'])=1./sqrt(1./ve(:,1).^2+1./ve(:,2).^2)*sqrt(2); - ve(:,1)=ve(:,1)/errfactor(1); - ve(:,2)=ve(:,2)/errfactor(2); + ve(:,1)=(ve(:,1)/errfactor(1)).^2; %use as weights the correct 1/sigma^2 + ve(:,2)=(ve(:,2)/errfactor(2)).^2; locs.(names{k})=sum(v./ve,2)./sum(1./ve,2); %average based on CRLB before weighting case 2 locs.(names{k})=v(:,1); From a6a9a985d02a278c3e5b0e8fa36f3312b903065e Mon Sep 17 00:00:00 2001 From: Jonas Ries Date: Wed, 16 Dec 2020 15:39:43 +0100 Subject: [PATCH 3/4] Cluster analysis added --- plugins/+Analyze/+cluster/ClusterStatistics.m | 93 +++++++++++++++++++ 1 file changed, 93 insertions(+) create mode 100644 plugins/+Analyze/+cluster/ClusterStatistics.m diff --git a/plugins/+Analyze/+cluster/ClusterStatistics.m b/plugins/+Analyze/+cluster/ClusterStatistics.m new file mode 100644 index 00000000..48799bfc --- /dev/null +++ b/plugins/+Analyze/+cluster/ClusterStatistics.m @@ -0,0 +1,93 @@ +classdef ClusterStatistics0)); + cind=1; + ax=obj.initaxis('boundaries'); + hold(ax,'off') + plot(ax,locs.xnm,locs.ynm,'b.'); + for k=1:length(clusterinds) + inc=locs.clusterindex==clusterinds(k); + if sum(inc) Date: Wed, 16 Dec 2020 19:01:10 +0100 Subject: [PATCH 4/4] intensity from Gauss --- .../+IntensityCalculator/roi2int_fitG.m | 163 ++++++++++-------- 1 file changed, 92 insertions(+), 71 deletions(-) diff --git a/plugins/+WorkflowModules/+IntensityCalculator/roi2int_fitG.m b/plugins/+WorkflowModules/+IntensityCalculator/roi2int_fitG.m index 8bd15aeb..83154314 100644 --- a/plugins/+WorkflowModules/+IntensityCalculator/roi2int_fitG.m +++ b/plugins/+WorkflowModules/+IntensityCalculator/roi2int_fitG.m @@ -62,79 +62,96 @@ function prerun(obj,p) PSFypix=PSFxpix; end nrange=-dn:dn; -if ~p.fitonbg%nargin<7||isempty(bgroi) + +switch p.fitmode.selection + + case 'fit' + if ~p.fitonbg%nargin<7||isempty(bgroi) - for k=1:sim(3) - gauss=make2DGaussfast(dx(k),dy(k),PSFxpix(k),PSFypix(k),nrange); -% gauss=exp((-(dx(k)-X).^2)/2/PSFxpix(k)^2-((dy(k)-Y).^2)/2/PSFypix(k)^2)/pi/PSFxpix(k)/PSFypix(k)/2; -% weights=sqrt(gauss); - Xmat=horzcat(gauss(:), gauss(:)*0+1); - roih=roi(mp(1)-dn:mp(1)+dn,mp(2)-dn:mp(2)+dn,k); - outp(k,1:2)=Xmat\roih(:); - - imgfit=gauss*outp(k,1)+outp(k,2); - resim=imgfit-roih; - res=sqrt(resim.^2./max(imgfit,1)); - outp(k,3)=mean(res(:)); - if 0% info.phot(k)>2500 %outp(k,1)>300 %any(roih(:))%p(k,1)>2500~ - outp(k,:) - outp(k,1)./info.phot(k) - figure(67) - subplot(2,2,1) - imagesc(-dn:dn,-dn:dn,roih); - hold on - plot(dx(k),dy(k),'+') - hold off - subplot(2,2,2); - imagesc(-dn:dn,-dn:dn,gauss*outp(k,1)+outp(k,2)) - hold on - plot(dx(k),dy(k),'+') - hold off - subplot(2,2,3); - imagesc(-dn:dn,-dn:dn,gauss*outp(k,1)+outp(k,2)-roih) - subplot(2,2,4); - imagesc(-dn:dn,-dn:dn,res) - waitforbuttonpress + for k=1:sim(3) + gauss=make2DGaussfast(dx(k),dy(k),PSFxpix(k),PSFypix(k),nrange); + % gauss=exp((-(dx(k)-X).^2)/2/PSFxpix(k)^2-((dy(k)-Y).^2)/2/PSFypix(k)^2)/pi/PSFxpix(k)/PSFypix(k)/2; + % weights=sqrt(gauss); + Xmat=horzcat(gauss(:), gauss(:)*0+1); + roih=roi(mp(1)-dn:mp(1)+dn,mp(2)-dn:mp(2)+dn,k); + outp(k,1:2)=Xmat\roih(:); + + imgfit=gauss*outp(k,1)+outp(k,2); + resim=imgfit-roih; + res=sqrt(resim.^2./max(imgfit,1)); + outp(k,3)=mean(res(:)); + if 0% info.phot(k)>2500 %outp(k,1)>300 %any(roih(:))%p(k,1)>2500~ + outp(k,:) + outp(k,1)./info.phot(k) + figure(67) + subplot(2,2,1) + imagesc(-dn:dn,-dn:dn,roih); + hold on + plot(dx(k),dy(k),'+') + hold off + subplot(2,2,2); + imagesc(-dn:dn,-dn:dn,gauss*outp(k,1)+outp(k,2)) + hold on + plot(dx(k),dy(k),'+') + hold off + subplot(2,2,3); + imagesc(-dn:dn,-dn:dn,gauss*outp(k,1)+outp(k,2)-roih) + subplot(2,2,4); + imagesc(-dn:dn,-dn:dn,res) + waitforbuttonpress + end + end + else %fit bg + % bgnorm=(2*dn+1)^2; + nrange=-dn:dn; + for k=1:sim(3) + bgh=bg(k); + % exponent=(-(dx(k)-X).^2)/2/PSFxpix(k)^2-((dy(k)-Y).^2)/2/PSFypix(k)^2; + % gauss=exp(exponent)/pi/PSFxpix(k)/PSFypix(k)/2; + gauss=make2DGaussfast(dx(k),dy(k),PSFxpix(k),PSFypix(k),nrange); + % weights=sqrt(gauss); + Xmat=horzcat(gauss(:)); + % bgh=bg(mp(1)-dn:mp(1)+dn,mp(2)-dn:mp(2)+dn,k); + roih=roi(mp(1)-dn:mp(1)+dn,mp(2)-dn:mp(2)+dn,k)-bgh; + outp(k,1)=Xmat\roih(:); + outp(k,2)=bgh; + imgfit=gauss*outp(k,1)+outp(k,2); + resim=imgfit-roih; + res=resim.^2./max(imgfit,1); + outp(k,3)=mean(res(:)); + % p(k,2)=(sum(bgh(:)))/bgnorm; + if 0%p(k,1)>2500~ + outp(k,:) + figure(67) + subplot(2,2,1) + imagesc(-dn:dn,-dn:dn,roih); + hold on + plot(dx(k),dy(k),'+') + hold off + subplot(2,2,2); + imagesc(-dn:dn,-dn:dn,gauss*outp(k,1)+outp(k,2)) + hold on + plot(dx(k),dy(k),'+') + hold off + subplot(2,2,3); + imagesc(-dn:dn,-dn:dn,gauss*outp(k,1)+outp(k,2)-roih) + waitforbuttonpress + end + end end - end -else %fit bg -% bgnorm=(2*dn+1)^2; - nrange=-dn:dn; - for k=1:sim(3) - bgh=bg(k); -% exponent=(-(dx(k)-X).^2)/2/PSFxpix(k)^2-((dy(k)-Y).^2)/2/PSFypix(k)^2; -% gauss=exp(exponent)/pi/PSFxpix(k)/PSFypix(k)/2; - gauss=make2DGaussfast(dx(k),dy(k),PSFxpix(k),PSFypix(k),nrange); -% weights=sqrt(gauss); - Xmat=horzcat(gauss(:)); -% bgh=bg(mp(1)-dn:mp(1)+dn,mp(2)-dn:mp(2)+dn,k); - roih=roi(mp(1)-dn:mp(1)+dn,mp(2)-dn:mp(2)+dn,k)-bgh; - outp(k,1)=Xmat\roih(:); - outp(k,2)=bgh; - imgfit=gauss*outp(k,1)+outp(k,2); - resim=imgfit-roih; - res=resim.^2./max(imgfit,1); - outp(k,3)=mean(res(:)); -% p(k,2)=(sum(bgh(:)))/bgnorm; - if 0%p(k,1)>2500~ - outp(k,:) - figure(67) - subplot(2,2,1) - imagesc(-dn:dn,-dn:dn,roih); - hold on - plot(dx(k),dy(k),'+') - hold off - subplot(2,2,2); - imagesc(-dn:dn,-dn:dn,gauss*outp(k,1)+outp(k,2)) - hold on - plot(dx(k),dy(k),'+') - hold off - subplot(2,2,3); - imagesc(-dn:dn,-dn:dn,gauss*outp(k,1)+outp(k,2)-roih) - waitforbuttonpress + case 'multiply' + for k=1:sim(3) + if p.fitonbg + bgh=bg(k); + else + bgh=0; + end + gauss=make2DGaussfast(dx(k),dy(k),PSFxpix(k),PSFypix(k),nrange); + roih=roi(mp(1)-dn:mp(1)+dn,mp(2)-dn:mp(2)+dn,k)-bgh; + outp(k,1)=sum(roih(:).*gauss(:))/sum(gauss(:)); + outp(k,2)=bgh; + outp(k,3)=0; end - end - end end @@ -158,9 +175,13 @@ function prerun(obj,p) pard.fitonbg.object=struct('Style','checkbox','String','subtract BG before fitting','Value',0); pard.fitonbg.position=[3,1]; -pard.fitonbg.Width=4; +pard.fitonbg.Width=3; pard.fitonbg.TooltipString='If selected, the background is subtracted and the fit is performed with an offset=0. Otherwise the background is a fit parameter.'; +pard.fitmode.object=struct('Style','popupmenu','String',{{'fit','multiply'}},'Value',1); +pard.fitmode.position=[3,4]; +pard.fitmode.Width=1; + info.prefix='fit'; info.name='fit'; info.fields={'fit_n','fit_bg','fit_res'};