From 31c577fc4739ce18ef453719ccbd25f988e306c9 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 30 Dec 2020 14:51:15 +0100 Subject: [PATCH] Add region growing on circles for point_set_2 --- .../examples/Shape_detection/CMakeLists.txt | 2 + .../examples/Shape_detection/data/circles.ply | Bin 0 -> 53064 bytes .../region_growing_circles_on_point_set_2.cpp | 104 ++++++ .../Region_growing_on_point_set.h | 1 + .../Least_squares_circle_fit_region.h | 340 ++++++++++++++++++ .../Least_squares_sphere_fit_region.h | 4 + 6 files changed, 451 insertions(+) create mode 100644 Shape_detection/examples/Shape_detection/data/circles.ply create mode 100644 Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp create mode 100644 Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_region.h diff --git a/Shape_detection/examples/Shape_detection/CMakeLists.txt b/Shape_detection/examples/Shape_detection/CMakeLists.txt index ec2042109ec..89925216167 100644 --- a/Shape_detection/examples/Shape_detection/CMakeLists.txt +++ b/Shape_detection/examples/Shape_detection/CMakeLists.txt @@ -23,6 +23,7 @@ if(TARGET CGAL::Eigen_support) create_single_source_cgal_program("region_growing_on_point_set_2.cpp") create_single_source_cgal_program("region_growing_on_point_set_3.cpp") create_single_source_cgal_program("region_growing_spheres_on_point_set_3.cpp") + create_single_source_cgal_program("region_growing_circles_on_point_set_2.cpp") create_single_source_cgal_program("region_growing_on_polygon_mesh.cpp") create_single_source_cgal_program("region_growing_with_custom_classes.cpp") @@ -36,6 +37,7 @@ if(TARGET CGAL::Eigen_support) region_growing_on_point_set_2 region_growing_on_point_set_3 region_growing_spheres_on_point_set_3 + region_growing_circles_on_point_set_2 region_growing_on_polygon_mesh region_growing_with_custom_classes shape_detection_basic_deprecated) diff --git a/Shape_detection/examples/Shape_detection/data/circles.ply b/Shape_detection/examples/Shape_detection/data/circles.ply new file mode 100644 index 0000000000000000000000000000000000000000..e876787906648af57071d4a3e59240753030d613 GIT binary patch literal 53064 zcmZU5cR1JI`@fQw-jO5{2@%O&8J&_!$jnHxG7FWAqLMvYM#(5MiR`_Pz4zYZHEdp_ zlKp#Kzw7IHx~@a`wKBG46XF*9Ip{1FzshQou*k9JCTc!Hp)(3X3O0x+_@qhmkf6ES7vkgGLm#i$R;`3k!=Shh4 zydlyBb)q`l`lpfoTg$62`rIdhS}IMvU=9(p*H=HXbkm~~es_Py#$!2&`R6(k;vVTf z2O^G4weZpoWA?_?j4$6^=!e;*z3#PMl%BK!bV|?zB)CCe2CxJ~3S6{z0 z5yNlq8BOtZ>4Q3bK9D4j#yf6OsXpu30@?4$R2Ls5Vg5zU4z+sM4Z@jUY|oa9r!YU; z|D`0Ki6z0u@tu=X>XjJYbFW(tF!sVD$AZ4Ohs|KPhZkF2WmXnotDRoTnj~TP@3GjY zwhTRxZnssE&?{u0ZfDY-9@Ghpd2YGBbXvgd_XK#y`!o+i2OA>`?K6l^&)elwzBMFx zTJfkhVg82=5*|%NdI|`%CiH*lK?ku;JB)W z^mjB_YTCX(0ROh8mA6x&@m8mEy5GtsLj8jhxp5w7-kXw{{kX(Bp@Ffm<3Iz#uQ#mQ zuM{qTKkm$&TK(-9AE)<+StF$fpu}M46XP5--oXtuSYV-SjoW+!U}BYlMCdL!b=^LkTA zDnFng3$yQDo<8xPd?)OYty;-BiTpP&+o(sVEP!_r;S$gMyD^o!2*0Nsu}!n1}eQ9OEBWHQ$pw*#lE9T5X*8g!Bv#{dDL}Nh!EfP3Kckh2m{Y z@d+)?vl}voWv(*ZL-rPC57R#N4uZ3MDtCn1<}f~HckT{pPIbYB;DL~WaHJ>veX&r{ z00-Jijs8-I&R~4P8H}frOS@rlfM8~;A+qPP@}EBVz7MEf{hd1W2F2UK-LBeqs|VpN zFD2i?8Z_Sjx(knTzb3&T@jrh{{HrlOp2yS^cx^jjWYFNP4t1m_3wMaL_Q`obuPpi1 zG?2?Gm2^KjH}>y%Q*d;Rsr6P>uiURc|IhsDQY94FrVIkfq>-~a189DQ?xycf>+OXp8M`}A z|1)0F;@eP{<~*pIuYSp62#xo{$JXKlQRMhPCF8 z4Rqk}Y#M~u9dLidP9vV;mw1-WLK3{n6-S7sE5pVsx9)exeS83NYxWe+96;;6e)}{( zH{CSEn?9#$#NjdfD^~Xl+BCbMqJI~d{e$?}&N@~og-ih4z)_vKB&74pq8up^a()Y~ z8USH85uc7~u7v|DBp7i1tlyEP1kAsrLAC34RzGYljhB3?i1fF5O0@QL8-uChZrmsD z&tl_EOnV=GV5Ap@x6&{@`bSTCihD}zpa6#Z=C;xHAw3m;o;sLAuD{B!|ENa&Q-Aje zya{e_odlb5?1x>buv^<<0@(wah!}azhwyD(zwkIPehvskrm`Vko`)!A{5TF)fk2<*WjhGWB5*po(@1~Ydi&qu z4PEdv?d@j`#}NOr`Hv5WepCYg+qYh(Mk4-0_k@iZ8T(;Cr-%LNg#&XZ%$#!^`vY%$20nP7qZJD)$7QKKe$^1h~9hzSv`)Zf_(-v6pStf#AvIpaX z%h*_$)f@p*A@T(s=TW|LjH6>ZdaVVD9{%~sffv~q?%Cb7kUj$3woX`v(<6KNm1w}D zQ3A(IJVesYx?%o3soR5=OXdMXl|OQiMR&K1K^mJZ`ID50 zfAGm0XN>c4Fqz+id1)EdwfK%4*}F<*Fkq8qcq=*><1>;rH)(jwL6?F?^l-ruv<{ zPRMvy&ico0G`}~(UfaloHv;XVpi^?>INFJ;VunH2GMRqpI`DptO9bW5PakfTnt($1 zS2Tps*@g0F-As(%WyVpc-L3WcvKNZC7ve4(=Db7jbDh^6>SB~1LLYc6_aAD8HIWjZ zGun{;<%&I5u0I7S<@(PS(RV= z>b}Cmu}KZNqs6}1bH6q7ih3o@!R0jeVCTn(kAiAj+~T!Dc$1szA^4hr`KMeT@)bTq z0@;U$o|c(3V}3Lo?(oJk5TI3J{dn6s^jyKyiQVm83-CPO?N~z&v|jjZ*c>m$=fglH z+Ffi3Wte|!o71&bL=sqJR0%SCnvL+{=hAN#2r%31mUl>UC59XL+)!J1LxN(po!@7S z(0(EAS!ar}Fag%hgHyD22)`cn;Lo5G395fgA=phfVtxir^tnvycfzrvMaJSj6#vG# zXSU46S^&+r7#65#AWL6n_6#Qa1AS!><-MS|kW(<8+?u?4gHxn<$>SXxqkbESJK&_@`34f|@WtK|DsbUpMD~ zZ0f^!6~;jfe;}{?q|Ub=%K5%If(u6LtCuf#fi9*LK6bz8%UzB5Un|s=yS*DiHDDw* zMgJV^#LrQ8@3qW7bD%Z*BhODU6weIb50>7hCBRsF)zQ)AXzaP~x0Bz+mn^`;9V)Q9 z2JOe&yp{GmyMc#`15}fq8ez!)In6ulW3w?}N*?C_ z^_PhfZW0N6e(Ts6Y?+SPzbqm0x(EzFWvWU_1v4~WF84?VLDnj0OjUnbARLW1nMSbZ z<;`(8@@~!`M+U7g^SSr;KfWVC%TirFT5q&32sL^wKQBjs)HqGvciO2Kp8){%1fG&$ z$Mc1eL*lgen%a>QgUfWAiu?D zgP-dtUulgbtL4C)hvbhjcAL|Zjer6FVfJ}AY>&E^=@$Y1N{JH*D5-{QXS^D-k z!mZE8<)qQ$q2;GXuf)Vqe*WtEvGp;_EIeKiS);l&hVgftBQ(9RAV5zS#;K;ZVr;zL zd+N=pw@L8sHLjW~*WY9KU8&mBGwfxsFDtJ1fOH_{$7;RSF2HyWJmG(o`7{8X2OjTL zN^ALA4#m1mW&`C=eM~nNpoN<=VEnkP@N>a&%)dAN^OKC7V{nR|N|MLA81sK_Yj~2q z5eJW$DGD1-pnl%MoHkUs91p!SR8=Kt8Zkfq2ka}3k^78mqlUyis%YIs)5M-iNhpBv zS*NaeKTE;vx1{%bMA?#nR)+e}z^*#%xpgL-r%%(g!p!z{{bEBjUr&^Ny6{{a2JK3w zhXx!`ed#j2Vl-`=4|7>#p$ZA<`K^*dt$m9G3=1k{9L+}__WcdMVr>}8Frg~e^hp{=^w5!rxz47WLFI05z%Au#l@o?@NC@CjYUJX>P| z)TN{moNa&ho4V8-JFN6HEMr1HPAdN-=JP)c$^O-a4Py! zT!onZxCs5lHgez3Akgb``}-%1&yui3Ys1hOytsOOEp#5`6VA$a#ghv-h+D7uCEtYV z^@C66uH$}}!CISv&Ls0z%>Oyp=BEduW`j8@5;Is)4_#$Bw3n=VIeZIvwm;C^i8^W}7U74oqVB z_o){*DJcu#oGz`#t@T9gxmM+Q!5;D?z*Jl93!BN{85!9c*vPrCX zp^1l~A)T`^M^T&?x1VY_Ml}a<{ZYdX(%z8a}Q-EUooi`YMVH;|Gtd!!`JnVlLmrBAZI|!y=bl ze&dqqvO14vK5Tgtm`Zya#m|M1F%2qmJV%VaU4*Vl7@w*)ls4}w$02KP-PY_$w7$5U zEUT&HCZKJu$#BLdny*78P10EwI5_bu+a-Xi1M{;ha74)_Y6@;yY)PGeM8NQ$*G69y z1$IJOd6CCCOGy7jmPh&3X-#0Te!wqM0L9P1*w?_5StD>_^Y^>UIjD}j623@xCuaak zz1chco_-$l-%B;i??l%E1H$AEfRo7H(#AxInRXO3Eqs1tx`y)Segp5KN3XkJlh4=M zU@a8Seq8~YO)si}Q~3{_?AM6T{XUm-1^aOD{^#{#+7UF~)Q@u*kUcNOXgRuQwkRVd~@ee1`$4Pcp?3Y|U zxc_R`(QyN$PY3H0|LWu>_+jIr^Oz#U{2XWNFB)kd11ZlhIWS$A#c+o*)|RO&1X$VS z=cd6>h4IOAt7B+&Awk8zZ#4)LsO|_)5d}U)bip=>82deU5FZfY?Ub~)8aTfSan27w zd=!4KFV=6D!dD^JTW&|8{awRjBGdTF95B)ty)-L`_IEvJEjKR95Fk}q$FU&3WQ@;= z+C7dTh9nqoF>!22v>xNXn#^cp<}(hzZ(g0i52N!gt}A*A?50GBdw*tES`*5*mhRnO zJ}lv&@wu_=fnZcOwLK>eaiL7^BuCi`V8I)tPs^)6(ogP>!Nbw< zyYUI7*tnjvfBD9LV;p8KWTgC>NAY&1&aTu#yAg`_9K69LkNOj!526z~~u;2Mf{0B^*%KKs$p+AXmb&NK3={YC5CD6$4Q(R zFpT&c-EGo&&07VHwnn7;9MiFReDm$vs8-Sx7}bg_*>#VA;g`?%9l}{RLEfW(CeN=U zJ+t()U)LTQ2g*!25_jlP{**WNXFD|B2Xkj*XyT<&eh|OWEadbj8@h?@=6^vz`C|;AjHD77wJFxsOM&p zQ6qR2D`@pV4(b0>q|lMtvjVc+?)c?!^*c7M3r~t9rvzqzteA`SX;U1AH<^welc$By zJ+xvxA96AYLKZ*pR{cf0ms z{&icfx*6XofDP{sCnQ^@WA+1#E8BdK>L8yy6SGSU8EG z#-iFRifRgbu1r0RlOj2vcQsjMJ^3e}GoDXp6<9czTQ8tBg&)uDOO zyf_kh{oEYz=J+{t_b$r+Z~mMLnJ2%0a<^SSE5C^R2We4qvrOWl<^^k+g2SkOm3g$W zR@`WZ$DgXs{heyYcz(Qeu3OoJ0Hi6?jAzMtVy7QlT64TmQauK zxm(HiI-b}~MIA=@ty6j?Q|)0b6raE6>6es>`AK>Py#ge#r5;QM)lo*v!Y=1 z<~$U6F7bjY7v=M7C0Wed_ixgxwC(dyuKL#Kd6`}e&(rkL7a`a6&>nY<3d&D#XWi9Y)Q%SjA>W6LxgTcL ziuSE5oo?{(&VbX&kO?{e|vl%{L}AAH1Zg^og%`i4>a5Nt~+AS z4R`17zx#p&Q%^THc~{k7&sFmc<4X}CKzm!~>#31wzv;Q~vk z>zFv91O9YouCQT2>;1dpO$DL<#)14jzc|K!`X5nCHNoMvF^CJXJjp+T?2jpS|LE`? zg4!CgU3p(ozPd3YcHlz=c|S6GF@lB%$#t;j zq2cvd*(MA>_n+lQUG)hNs?a6Tv4HC7*4<$-!RBEYIWy04=^~2f<*gfcG+k@pN#Yro zXRlB^_n%TTm3ub;JE!}LqQa2=T$eQkFD=Ia|JljPKWk{c(1~-u;xg!h>P|fZhrCf- z$*TEN>a7fdw7)WAu61a=kGiy5XPhO%h^Or2JR4;V|P@Jm@oX9+r-wzyn6|^jGqWUHEGWWNdXgBPA zG0Fe&8j7p6vUg#q*u# z4UqmPIj%1`MGZk+`2wjUQAmI0^_%qS4?CceT|Y%#Oe^MR`L~nkSqdCb`&xBjj{~X~ zZ&_L$zE#)4P_0mt$|*EovI3UpzkVhGYJSd)!((W^Bx|)Q{JomszN_D^DD&51{v}jJ z#b{Gz0AF&;WNPv_hP%J-=2t5ofE0V(lYN}gc#kl-GdPew0cXv{J(d65zpDDrxRBDM z5gt^vOn9wPjrq}`Ad2iAngz$Z=uBTfAH{I@`;^>*cD1m*?|a;3pEL}2xUl)eXqg1O zxE4BttdcQYd&Z^7nS8&;V2MB3PQ3!|=zM6y@`kDz9%RNUIOaW>#_*M^#FM03iPcG`9j? zlJlzeZv6>T&;n*J`&;mS(4>TU_hI`xWoA7~Xz0bK`Dcs8G2A2<+ zw@RNt`-{KbsTa>~bwYt#W2e|5TEE;qn`ceDhk-%$QCmS%w0;Sz)mb-%2O;f~w=rHW zD4tc!pG;|tIm4uAzkv?&J**x56ZFdJmEX5PuKPcYUPd=zc%5D2mtGjp86WaPr~0@ z@)Go;kZX6x?Mth(7=KsIlr=?i-R*rbH~!QR^VP^^uN1GoMEQLGc=ZdC7!EkqH?E&lL;3uX@n2vyJm8s<#+f( z{k=^Zc^|eDC%ZC~ZfCYM!3KM}z|wG(uY&Q9I1U=m0JD7$s7zRJn0=qjVBIOBeBV(zn6zW~yU82p9+2;cGaP(QRIfq) zSJRaqeN!HWHh%I7Y){boy;+F=kJ5My#${@~TXR75{lwwMqs4x`FxW}@d~Z3b?`e+; zHyY#WV19SvecOI?KJ)JDpKEnF3xGXcE)}hY<6y(01h5u+Jguom^R+uP^onA5H>@I% zJW4(x{yN@UNy6Pdpx0~d#5@7a2UFbsnf_wEahU4c*gt$jd!Dx|{sn1Ypx5Z6iKr zX;p1ZE?scDd5~I69p!&D(YgA^bVOK_c!4(UP65WpC)edK-E|T?Vvu%)n;OMcaKc?` zL-P5#t>nw2G7P96P+|C%K)*c;xemYO7urVWA%}PAc;%XQK$T0M2Hr_kV*GcVx(tLF zaDeZY((%7$GZ?x!kw9fcSF?MX=VF_QATJ*B6QINKbBu;L*(mJp6k@ z|68a5S}(N>15f;4&%nU~4)L)ORNpI9g=w8$k3hwoZte%i(fo#GoVg>L(gcgg{97vj zdH=}Hq|8*%P!F5*ysGtPzF|Cnbgmp#^q2>hk(Re+i_tk1)#Nk9GnsfOyvFIznO%qJ zdFs%V*XMWBuzmc8T>E7_h6~Xfeok8+h6DU_pS~Qc!~EwR%E&i?wQx=A9jJXw!f?R` z9`5by_3%04!ruJc7z`h{_`rvgdI5aor#!Z~zX!v^Z*?ehX|}@q5@)NVC(t>DN4aLZ zYULy#kqF;j+oFD{HgmbFi`;L7+41e0Tu1~v>-{mluJ|MQJ^G98?;Gu3JTZLtpS@Cy z)I=!Rll4$+#Tk3<-Mp`Fyv~xKA4iTAH*YhBGr1d2ST)>nNnDS#%d?LFFG9G+h`K4H-Y?Nf~3IsHNJ#|zqXL+)Le0FCy9 zhX?mgV*KeVSaeQLj6*$F+BfHeQU14jHl|_1fP;INy|{YKQ2q}O(oZgSA+O)lzP#^L z%Q5~G`{sNMyXV1GMbXQB_zX5jP?+`)?|7YC7caCEtbgUSaFc!uX4q&(=3l)x&U!zK0P9(f*fh zwcdlIx&R`-D%C`}cVT|+J}Mnq^%#fzp11W1uA}+to^tm!AnR=MOO$PwHJUHmDYn?7 zZwPSg&2U%VW(DRa-kp=dn{ggqe$A|X^y>hI_a8XR{LsD*PE1oDe5;MlRYv99(Y5r;0AMLq(H^Hz)^C|G78prQ=}-0{C{D@WwC#^`Gi*C066f=V|hP zG%Qyw(RzR3a>|)Fln*T3dj4GPNBMk(%jD`-O&eT#&UfiXD(aW~p7^LU*b#xzg8LVi zG87N{R7xMN@)MzGkr#>ML=ZNvJSIA=A0;HXD_`!7(`Dpmzf1x51Lit->oJq4ybH>= zS8RCn-LH_({~f{$83maIbOpwT+{QKc5nvjM~{|C^;m4GUhDepPbHm zuAE#&OmF9fp`Ks@t%IHV=e1DYA^HcX-ixA+3gOAG6wLl@#Uze^9@BWzpX-fir=sP3*?LBi%uF=#(1&)aB z5SvK2lTThhsIF(Ywg4<^!Zn?XyD(g9Wmn)4u`x)oJFR}&7p>oBeu*#UtP@4KE8@bL4<;-caS%#YK2 zChvQ4f39Bs_-4R9w4bGO|Lf;chlh@Hb~i{|)0n-ma9h#E`;Cy#%HqILW2FDm$n5*4 zowHz{;3lyo9PO`7svlY4>$~Ch)xQBp{E0iz5lG_{d0cB|GP26{bDyfIZ||EO9j>U zr4^O_T#Y`ktEKjZW*DmPkLCtqjTwp1I)5g-(+RzA68kbzUq-%X7c7#^J8`o~0e z3|_+>+6~>&eC>B~n|N#71K%6rI4aiWF?;v@k$W%255tFB+{&Rss2>=(K-#S+QVk6c z>>C+)M}j-!^$e7yRpuXsZG$p3Y75B!4|a>AFW1}PZS$!#E%JBRJNCbiq>I>FG{EgK z&IJC)l^CutUSVuCHU~=a&%J0q3}SfPfl|D<75V!wxrK->Bjn%k^zKjd&Argw+(K18 z8}YwoMiFnbG6*wcin;m%5&zQ6Gln)R$*?u*pDry95uv%r{UYk=lHePME-(AdtL?|2^czH=G7{k|i)SsGE zkne|BcJQ@mqdah*gP6YWXfLEM6Moag+6>#S9o$>Z$vC6o=W% ze{g6Y6c_?}JWj=Q;_l=bEidYe$0}{)ww~6o=}Fpo?cdE96HihMP4j zeW~#$0GDZ&rb~XPzL@`bBmE;44>~n|P`7NM@h(=ZTvonN211_Mi!PAs@Q(jGzi+U` za!r7y!(&5whf#mqRhfM{?(rnJsy`-NT}#BCo4$P^JnJzL921{@n{_t>!=w51E#e5{ zKt*u>gYQ15-#hF)jQ8#u2etYfey>-Q$dH+&swMGGZ#UH@h3+xO4VAK7s5X;YVFrA}7GHyISfjv#>l^DIu}VGK}Wy#r5CZ`&7rlF0btH9xYVAoG(4; zojx}VCzU)7`?}0wd{!o7PmyHDz$8eE#hpieSL#^#R?VLYm_@RNXK$c-ec?2RsqJqf zaHMw6W6VVM38@1laS{<2HpJvWWz_G*+b6kyxIPA=-!;lkWF}$bk`tHwv**VId}~^C zNQC^m2s=8woGEG<(HH^LjouQw@1D_>T1|9^=QF^2f-S`K6A*v?8LW#jMG)%nuD?*bvV-V z>-}9EHotHnBN}&p-~`gYS>}o)w|5f6tt{-O_>Jswc8#_FPyD>!H@QEUVZ9AuuN*%OMdCTVH_f{DXtZZWV+{yR6cJz6Ek9NRE zWD@)tR}LA8`he+EcfGdDHG>3x7=>#+O+ovngp4KqcSJl$->p$uYl+Sky;W|#2>D$D zDBd0G`GxAjDD<aQ zP$AAw_DvYl|I&3!9gFlaK*14{!%$t0@#!(nT`7Nphbc2JOj0dTKTNAoR>2xQ4)}H} zPB+Y;_0GHyMdT9hgJQh34O>z3m>-+Ptv#JxIB;J?+sc3r@s}aVF2(K1f$cfDLpuN5 zfBT#uPa!Wf3T#K{r-)LaU}xSPGS7e9b6^f86?^ar431#rp&7w5zO_WdV(s z##FaU*dhe-#3&p*NIr+y@zdPN%*cFf9Q4$sUC+})`4BuDyfv&p2FIC8)DOLy#q3Q^ zO03)TjDXl72gPECT)5*$DcdqK-)0^jHRrM(V@LNfGPPXcY#KR!VnRE)|Ab?H+RL02 zhZAOj3`gc^aVE0f|M$-#cUZU}cNBQ86FC?1QNJN_oY8w}V;bh~bDCZ4V$<~J_fn7$qo@1h@j(a-p`Eu#pBoG&z=}}60gW1m*rPsSW zB7tAbF)?KWwV3{m%|TC3rB8seQxVqhs8KzQ*5B{(BykpGJnQ)D!Z3#UDa@dA2+1V? zeY2x0Ba=w~%(5|YK9xRTcJqm>+&{nP5QzI>lEY5`_NL%~+h>&j>*`;4rIK0!qkk&R z+y<%_P2GpG4W}mnz2#^B-Oo|Hl~62MGCUxkC$GE~6Ba=6W~}qEQv1v#(7t@`;5S<@ zjA!iu9-V{ad!Gi{$+-y@jo5gdD6Ou9DU##u;14z1y%8AhcKmCCU;HeL-`36(WkK;- zTCTzWgmMxzm<7u%8fRnvyEB;2<0KY zu%?}DoH^r$@z1`&T`t=-5Ao7=mwL~k^SYWP?TTj-lOXy^#Y&-O62|iv%}E=VHWE1J z*DmAfi}H3{$Nhq5u|yD?${M59hJTRYawU5i@7+tkJwr7MeN~c52TI&$!ko3%3V(}{!WbhW_1oNzz(4L_Tlax z4A)BP$h=)N2COA*WIbI_z0OW>8?0l&!-ymo&$hiUChA@1Ea#w*=$T(07kvL&KfYu9*_&SY8b07S8 z?#{;+SCp?lAFJ-mawdSYe=2$lFC#s_6ZnI@V!FWJ6eiXTMYLXajg`?%o}B>vD>iSbORq$1sqgt9C^27azlqY&=0zRT=U4F8S z`tvWl%;}bb6JX(jOzbQ>nwPn!!?OGXRWO1*;zK|X!e>*?$uYqS*f-(rkoOJwFRJK$ zMwLwd9SkdjQF$AbPeRgXVx4De!9Mn`tImn&xmOa;xY2)GAip1!>}@db!1RAD;&O^f zsRp|IC#}3xGKcBR`f#pQcyS79)ec&z{zt&@HPXIJGfCDZ_q8II$B9S z$IjPQp~j(jJ`)*P)-6^GYTlQxt)4-AmMq`cerlWt51zmKBB+b%sl?v?bCim;z#=SG zg!-R*wk$~?YzRLOz7D+?5VwH*T zS(?SiDH-3vOO{Qk(~Bq`b`81YP(^fr&?7yjLXilsI@u~-1^}AgO-47I>*zbknlhkEdQ*TCuhmNQ92{%q+@xwwd1RNCZE@T!|VLITh4%EMAE&xlMFr|uw#;a=G z@s59O3@EW?{M^Sji}|nIFLItTrVebs3Rb^Wiu|kf&$~P%-}@Y1G*K0KRfGB0{1EY) zBBd4>;xF}$lJD#7=oabnLU2KE0W7gh%ko`D@8f!D6SA*HSAvtrS>)?OQM@^}$7)Wy zOatL_rtyAyX#ZFg8 zYBC?d{EvL)bX_$rgXE)z(pg5duGM5qK5%!9L%ZEQ2VoG>K@^5kF=W_B6aZHMf6gmMmNR+K~h>_g~CK5y9J^RZj5=>hq7J*Y05ojGVSjq#~A z$$f0!7YB@gL{jOpq50Jw@htt%xE1g60l&fpfjXd+L!ko&~m--d@cnDjxH4aLtD+I_JVTuqQA`^{G#1C)RIjT2@MkhxR3n%A*Q?Sh%tI?m_5hCBCglp zYT<>?L!9Cxs9rOt>AhnOnt?O!q*8|Nam-#tFHMk5xCZv$FibQxLiWL%U7S1}vE%@0WITkV>6uGN-Quear;i+JEc~Y_Hz?(K!VsppV!mL1cfvK*N6VaW$}VB3=J2 zjpBJC%`9PmGWoovwnhIS-5BOSrHdjfU-T<1s9Aro#EaHxzkPPxaAO%{tn4|)TYzwy z_AV~r=c$k{@25Bq6N-l~;g=3D^EAW9YSG~1S+u^?zOo6d>`Dc()b5ugHIRLtjqqi* z=5EmVxFhsG@^?i$<8`*SYRlBghw=pJNTXyF4{4HgBWd3Ku*&7DlDIDF6OJk!%GrR~ zfI)f5KGPG~bDfDFa?2e8+HS2?N?~XnB|AhJRS4FA6UQch91lY4(Nt$s_I=|xZdDG zx`puXd@pECE!-8^t1><5t|^&=1$$0$-Qqy}SwgZt-tR4fGkfcGRnw6E1S_w*$BTzy zrJ{8Kk2YGr<$FITo|&zNzpVPTek!8)&k&jbU>x+M+efAasj z?sy*K-byg)b2VpAJPGsf*Lk#NbvYL(WHcCv8zMhbcKv&a`lDcl#rES80gcOuJ{Yr)xw44BuN2jWPh-8#U-q82Iib# zE))nw`Fzf{-swVqB}n@~bRGym_9y8y`}0GmL2c!qRs&HKZ{eRjJ3d~n1vvV3ZlCi= z&xZ~*2?hQ0;AtdD+?xjNk4&yF`~Cb-3|dYEo{pAAes~(|OphMHgXdAV4wFCe$Mnf{ z*5q{1tb~dC|6~6ih3xGDw9h^>AVO!O&Hk$PNzA^}C$qV7e?ClG_Bd+ofZ}aDrG%lB zyB{jjS-PcK&tvv~|1b_)ROP^tC3`}8Gt!^F`}ta;O%L2(D>(S_FIw-Hh=OZeKO4#4 zsX6S*wMO~BMUnnhbusxq>a3H^gS%y5Cm!m#Z@eMDAccOt$Ko}}zlXYmKYG@jvza>s z8*5*TBz{5qTZg{jpJ)z&Gar`EnJyr_YJjz)vhzEv>-f7t{tgx6Z#j`@smNIay{k4W zi1o<-g|y&h&+;jlG1q3=bq&SCp_fJfvFkQM1@6P0&&5$)+oz*8=uSy~e|?+jVy15^ z#7o(x*c7c}4)w2t0kFf`b17dM+_{lhgQaEz}_Pnz>h^=c+xR zGGdXMffHWw%kTgBz1dV8ed4F5bx`$4(g&AcD9)2TmnCc7&BCYV-uv zhOZ#kXLVTPG@4((nbx%{KU#tEz6}A@T}Yp|ksB6o4%7qQxtHe|AEJJMlkxvCb=~n? zf8U?7LZqyuNGYpj6_u+*Dx{(kAv+l*D=m8~o3dwR6Dr}_d+)vH+uQrihY+dn>-Tv4 zysr13=W!nA{k-Qs@44rW*STLwBG`V8&6{8!rN|K6pIoe&KJFR?us_|@w)-7w@9v4< z--QeP@W%7~<6L#9y{!F)rnM5KAY-YjIs6*p$xe|6?7xGifU(>CFBj6$xG3?mjoF+l zfrqn_ULK=F_VsEOErlG$;B(sO3|(H-FNv=R*Y05J4PE#B7H{|v|4%NScPe;11}{H4 zKe{GN!0`|2c}KeP)qshg|I(&<(f{tc>hnTzh6KR*N1g4#D1S(;{~b5tk`0Svd$mV| zkw587ZcyA)ZH3m~PT%lqLi6Kbwur3`6IOqiHi9d~1>yg*`aAbZ-2mhmb{f7DjQlg8 zA6c%=odL`Z_XMmk}p}t|pv63qtAN-$rM_ z9Lld5>_<_1FZ0WaTV7~^?M(ggl*U)8FP zd35my?tiP##-)D3IUp|exRgmM^3N`?yO=w&7f4R6p76PYc(RtbD6{tB9x&zZB~WZz zzo^{{8Q!SQ0)+Nq(?=vUU#E0j_;%+_BsiWR&e^;l$q&YktkWgcfYFDarhgbAxq>Pi z_fT;?tfUU;oP;P|P2?Mfvb`lkees+mz33|3KQ+%DT}=^A06IlnzRBD6PgMV06h0f! z1gQO8j5hY7@wG!|MX@H+4Y(ezUmF^}su(7^`xN?9+MPiHfeShllSClb(G=_z!S4`bBz?;lV+NqXik2 zIR16>f}zS2&tRGjSETVd^uM<*v;AlBH31j}cFmEtov(GjI@2&b^$jX<8eNEfj_@4s z(2!mK+XD@CHAA#)(0s9L5g$k^g7%4nHUgKs7$f%3u{Fa2`!w!{w<#lf+!Bgt5QPH_*eoX!x{G<85a;Rff zuKVT@;^AD4ko_@iJ;E_NXC=vtcs^<7e8-Jg1C-*W+xCki|9_QM@xF|G$90yA&3x$W zcijJ4DN7^twV&W+A&S@9+xCb5?5{D}joqK{++X4_`yTa6dn#S{HCp7XlBd4)C;3tF!)+8g){4%?4B&wb<3kc~Ve3l&JxIbwO zQ-=CN>tP4Y&;NdjqIT&NW=1GU#s^uJT!XaR6{ z*P+}=jrxn=P14o=FbZDTnn|_&AmH}-xKH@6s3w3r554wv&>}lV!ft0-i<^M^;V=AW z_aHgn+I@h|75ted`R$(=CR<9Z-;RG-G*l*ZErWg1t-U}|HRSw zm3CC6@_oMm_dE_}X>>;CL1DTp$0Mi4VJ6$ZOde*0v&tvY;<(lvd@h%|U~~uBkKgev zEDC$SCZ@|il%g8NXR}+%R3A$h;d1|gyOA-<+i$&EzL{B1hIfVic4mwf;y4+M`_*5( zpN5W->4}PKRk-}&r+nEM?0i}9`sK3I>?q%$th?%<@qPx%pW0* z1D1cIeZQD7rl6#zDUjfLrstd$vLBRbeEx9dBJ{tm!zEmc&Ow?=35~TRGJNaGvVT`7 z7C-*~=XhLnlre-jzBq4y(rU9bHn9|3jHIwSArtU{}V^fmJRv zuf@w%YWNQKfwQlD-gnWXxWW3Xh`Q|3G>qpOj}B5p{NdQ`M7kg{4;Ir~#IGs$fh{~i zQGZgz^B162z@xcCuE?Jy#(yIU7O? zP|w$i*69K|cc#5)Il%pA820I8uz&o9@{l+t++UNJ-qI>=Td#{TW zX7#X{62+g+!nx(+Sbwp337?|lZN>4t-OHHING3z)mP+kB8+4AWNh;iJ|78~HIrzSF zZb18@ohms*8HrIK{bTxe{yK^e6P-VQy_lMT7Sfl*Ildr&uAM(0fA!x4@G?Eu*t6|j zwzYjHC@Xf(K$7R;bV3!vPa$hTu z(Ftipz{B}dm1Kg}*8{7ydtt>QG!&3H95IaeDG|Va@;7$A;5J`2?5c*&7kUoU2u_tR z!V!%Y!Qugwm!7!`|f$x*>w>vxdpni zY`gbvzAtLI`^h}C;8zO%B8>b-yH8T=$HQ6p^{XFSzai=`8)Rk6yfVDd4+v9vo$WnH3Ynp2FXBL?v0}FlW$;4F0kizo{Ib;{3BtT)WL{& z`|g`tJg2deELOs z@3Ya!MR+cID(7TrPx_0)G02=RE2S~X!?qs+=FTl#T zP>Hvz&=peD=T5=Abla=nrx8Cp5-$wDoh5@bOZVR2zYsrvl^FT;e4oSKuhjjN8;tH# z@t-RGue_@sSe;<4oH~!@-z=kD_qE0s;G^5r5nnn`JP+vhFLD-`hMzNHIBnTbJa3yH zTsuUu0Czrbcw1nN<_kT~olxx400cEf#KTc&zMy);Et9V`3*Eb2&5aBY-wYKKl{mhR z0ng)QT_vlCZ=Zz{^L5h}pnp~6%lTx)|3@Z!C~c`Hp_R4!TRUFFe~;cnHk^rz@G^Bk zbY?QTM^r;+5p*z_1YZMS{BGO%d+v%@a7WrK3u4Hp@vs+9H4V_m8c({R`&;OZCO)At}>F%I-lft{ zGqE45OUbZ9`1;dy#OK+MhnHmA%0SZi#m~#z&TB+6XHv315n*zK+UK>G$o@lX%gAY; z9N^}$dgcWI@qs8;QXR5?8eZ75^gP3%3&#^($aaQ(c^(k*iM1lz>g1a#&N7|Jo`Rjv zEYp5HLi`+6$-Zui<&i?yDJsKw<8V9l3x~5krl%pzO1ef>Zv`&DP!=y2urv$yT&tm+ z`aTA>{9lo3E_S?h0jBQN8|UXJ#r4%LaqblAn1pfOm+$R2LV5DVi>EFVvG_czR1!Hj zh~m|fsajSw7N6T&86@w2MDdw|BcxinY7vS~{jsexL-@s&53$>?lcCkH-H+tvc>I4m z7R`F5zAVD#*{Y~_>WFVb`I4$*!z6f)BS3Sd4gGI>mZzX+;UY9qq<&eyf%3f~$tYc! z4`fKlu4Ju|Kz^_qle@yFzX(a4uWRncRpb9Fb=^U?hKC3#XB7{RdJN(IF%CJ>za&bA zJ3!38&%B5aVs85qGqLw?S_FcoYZcM_a(C0_sN=Uqc=B@B(uH45xP8gv|6M<2jLmmS z6aUg9YH<6)h7UeeWz0ZJmqnJYtH}Ru^zJ)!S*O7DaKjNFLFE7QT{Kq`Ij5m3W80y= zakNf({dD$oHMTGL*;u~ns96_o|AWhip18Pqc;P}>(Xcb(XUurQ0k)gD@KW=r#>pwf zPx(pT=~LLee1(JgUY0c)kB@uzRaMmYL;YhSoeb@0Jigy-xv6DJgk4@Ei&{mi=oO?4?>EDFnKbh5rW>{enHtB0Kyy8Xk*BI^5f)JWH81N*}B!~~qUk>s^5q{$Y zs7p_Ke$xl#dnKEvm^>u9K<4{TuNQWpd`}(RH=r<>h66WuFG!!Q!TtaI_rH`eCn8v` z>H5mBZxHv7r9<}Ia4eR;HCN7rZhK#;yzSpOi#rL1)YIhdJd1es-t_%|{Hb~9{bcq2 zKMmwhl{dGO(q&wsvPO~Thi&l$KIMx}q@0D;Y4oo>cuR0RBB^iq`0OWuQn`y>OE{YE z0t)<;cL`5J#i^L9VeKeBe2shWD385!#=&+w;sxlz?TC?QUvc^@!1aCn_5@ZG|F4}? zd2zd7941go9AFeh@jtohtW1;g0^~iIH~P37=`Y@LKkIgX6s}*}!NNX)^2zYjnCG&A zB2WF6utlfvS6Rpl+eVwl!u1*H?5%)QL!3fWw=f zKIbYN)PUM6X{mXUcs><=2?{-Ayo}o0y2#+OBZmxq!MvSvGHRFM zBmPrek-gB@o;6!<+y35=h`Mt7B@(2SGS#<>LiQzYlOIS#RKZonNLFcVA8yNUq_A&1 zT0f?sow4CTO=aZI4}1JttJ}z+KHhTT$pFfWc2aOuNHR>r9k)Gy@d+dTB)m8+AB=rB zRVTqmZ~MCgw;yFudsx199{%_4l41x%cwQHd#0A?IK>apaOZEJ!o zy*Lf&zqv1XPmiJ;o>O6tIlt{Z3tK*p>t;xXORUXBQF0D>Wvbtnrfz(EU%oV!iH9+>X!A z@0D{$M_?137TvD6Ik4qViy%-Vdu0)F=_qlf{6l#lMJ(${FCqbQy0a7;+(7ZiGApc{ z^3x(zujwKF!QTDbvg0n9zQS5gg73AQn_1$}I!?}S_l-OBv+!;K?ZYpVh@Xk6j~5n{ zCO}?_wX|~~;%Di_ub;1D>&>C1>Z=S~DBrt$D0^%c1j|39d{B$lYbi~(5r3Bn(Aq?xS5^bXPkB!7y1yQsAmH`! zGq3qj{On~a{n3y0m()=M%JK-*uI?a}hP_x_GqZb_CSPku;Ql;F8Izn%zX(qS$XeVy zg7RA%Rr#vE%o(^Z{*c{mtqB}Ig>SJ|sQL_~Ux?!z@Im^$G|qeXv`hnzSnnYn-AP=( z*!<8{pZa-Noc+)2z&3uq|LY1_+Y8{J)lv27HpEZczFW_1u=S?>-+R7?$cS&R_bG*Q zWBXmIf;aQddVR(1kL9aT2Yn&J$U_VQ;|_>#=h$6a_g=~YQ+ zlN*6f>(Qvc&WCR>+Aj5hTdX^Q6wN$t-}Z*vMjH0MPIBzSd#Bb=yrpCL$;stDg4K%) zS-)?L`uE~}4MKC+9E^P&c7DYQjjzNV2`T$cnt*=W-#kY))LvE)Ct^fNfLB@@_>M>* zJaH}KAL6k(G6!vgF8}(A>^Ds@pEd?FknWd5l!*rNe@zmzp_$AKa0#%sdsWE;f8*Gw5SGp;q#!NHrPsL>VAYMTdw9k@~ zSWH5ACi81piZd(Vkm7{+&0yr8u*gNaeLtt5MAfajr)DUAzToiCl&vF!wa>3Q))b0x z`?=ZYp1zNngEyqrB=c>N{ZpQ&O6M9I!P8oTd&ePUKj7Tv-bdJcJvtp;AlHfFZF8jh zpnt?XNVTh^ZoG=-wM%tse|PsSz;6_$?kpR~|Nl}qQs?TY;PBL7{H8Tp4;C4cJ9DCE zU@jGx!nzyE4`;q#uYX`L1wyqK_~i!Bxt=wfdE$9&{-*KbFkkmVYe90g$-VWgr>J*In**pc^ z7Ypvh55FaC+yAFi7MhZBnSo2wJV#&7BR&W<-p`|+p8{86coNRsL40VOX8-PTpbR7w z`M;`#C~th|{B&WEZ3;X(^&pQ2+gI4aIdO45*1|al$hIwi^tVUr2r0plVh{Q@fYlUP zuKtR~?_E0O<|wvom{D&j>lnR&_e&$?&+{%aqj1g6=Xc#Q0hgDYOrTG-DT5J%)DQQI zqkg$i)O=NIlnetzBm27#<>Po9_j<`eav~HvP<@9w4dHn%$?a$KqX!o4>8+x1Me#H5 z;LVgnr%T|oBW{#Ayr>`757|U>WA{b#sQZt){2jpYbS$1YlB5Tr@3Z*^H%+9kpJYY* zr+yCluq4pz&Kt({uYNtPOP^l_(jI@7D56CAvguMDe4|8AG-*6U8_xl&+Sbh8MRu{n;#Gg79rqstT%3yZo*?Xa!Xk61M8c~l= zVdrZh$EOo%b8&xiMUR-nix96qg!M}1PtF6axsTLrtthVEY7I`S*DQtR z)eZJ5mMD&$q&VTJgXPtxr>PDubff*r#d}`s30`F|`rwL0)V6yF54i=GI{%Vk)rF>V z8@T}7|Nc9m`N8&|P@6Mo@zoCGpQDo>Mji^!LhlD6KfBA&{FQLOF{t{_R}fpTKmEu9 z`Dfp|E5eiKI)JsM&BL2o$Uk>whQnVeX9ER=TABU7kp1L?=L$kEw*ZX+t|8GS)Go&R z2k%v%D}u+`&in^!5I?!kQ0iS1BS4pWgJ4x1w4cK$qrB6G$pg%O5u-2IwTS!KQTuq2 zCVMWF@C(~nyluTNv*=3-AIkwFRC|Yxze4!m_x*kHJRG~v_Ax{|BpTswv&h&ns+0$r z{G;=wUn3s2%!!yyJ(z$m&$uknWTX6cQ?>JaA7eIHh)xaSUqS8l{#`C@T-6M?4b3J4 z|DpDZY1c9Gyo?74c@KSEYbC_y{Dqie=cWJWW5SO9xS?Nc3*mQ(?_iA^8Dh0x+8s^bCGQzLu$-@Itt;xJZKWiUR3}ee>1#3Ac*2al|0pfK7$!p z#8>aNl!xMj+|uEX+}4$V(G1RB#J-={Y8TBL>dV1(18fO$a(+iyHs*lAMi ziuQ|oLS*!mvG1F#Jo>B%KahUpC$EVABIiNZ{xv72OFg*#XyG7n!)s+QB!qpc`7QFB z)x*VR>R2q_9%=gkyvlKX-GGaVALt*$Rx!@;uG^@;3XY^rpBkuulhp1yF9T8kM(giA z!YPccH}-YaJ=zxE!VU>rKUnqvXZPxOKgHU;<+nhQeEkya9V;3o^^4hs$e;V!AA2Vl z5P%Qc^`l2P(RgVV8{_ghSPj_d6oT5{qxfSf^F(j##XQiwG5k7e7Ugfe!a|llPSK$E zL4Ud@0~%kmvH`D!GGYMNnY@Fnf6)KBua69cIlBQGM@MQFdDPwk^BHwJ={z_*zVlM; zIKrR!(y!#4U@~+S3^F&)Mf#>UJY(`>d*SM-cXM91(fsH(q z^=~^vF$Jpf?0UalK>0xXj?7EEe&3V!graq|3?NYInH`XuhMb-w!@DTWbfxgWt!L5 zg>?`f(T`tM35OAgXkFt$?xuZQ8cw>A&k`{!^x7c^F97^$*g*fEK39AA)rN>a#r zlJy9TEvX9aeunh935-{Ruatp+1^)liYY|QjLd8S-+f!g~@6i5&dF1~N`GIdvb=x;R$U|bBjCV=`Yhh<3RNfimC(?UX?S1aCLm7h;;?>UkMNAe>Igy<0d z{dZ>UJh69lf-L4_dwBQ$r4FQjV7_L$uXO+>y?k$O z8i?!{ng+a?YKjEi|CoP;Z;Njgr@pP~d=3GTm8=s&Se~-wPfo6L-M^nEz-0LjD!F)M zzoup6)N*t;6xHdhi7!I&p~`SzzdA($RK56ESB4+?N9coe(0to8ERl-2d-5UjPoY?j zR2lX@+AiZ46At?j-*zW9{yqQSG`!C3rNSbO_Gh{|ZbZ|ctpa5AVg^B3)Q{}dPC0q} z*!K&c^m!?)hH(B^TJ~Ps-SGo7^|c%k?L_#!Sg&kEbq@d+esFcgANgl_`Ec+~)?g@c z-7j`z3E@{z=@RbCXn}Tde!r_H5O3YrZMp`r`_e>$4fTCnWWQsCc>GQlmJh65dNNUh z@)-}_%)g4*_k9hmwW-`#zQ5(?gL1rd6xe+{Kex|`Az=x)+?%%G$<3iGP|rln9^01B zXuC$gq`lY-CWl17-k?K%7@A%cX}OvVOrM*H?YxHY*re73g?=dn5%i*~uNzSRPJl;& zcZPpLBKVuMsfPSdzf1A6uJ|m>h-Q_#|7{%SLp*Jr$eWrXNTw_4bK6#DfU98rj(QXU zT1s7LQLh-q^(Q0#Q+aRm9i}(^y8aL#f3B*@bH`?4`^@n&`GEmwKOw8yB&9qk589QH zyu)45cz=^>+A02i0-l#V=i%*!^1Zj=U4FHK)gV#zt`O#D+;6+g8^xA==0TOBz24&> zlt1JL1g?Z$DFHSPk(y65k^R?4|2!Dc9tU15N>WkwvpD`h5&zw5i&a3C&n39X7R3kF z3EGMvsd?~{Ly4=71+AlVOVe|2RJ(zuI%c``E6D$zfA(Chq|Jl$ZxbVKZmVz8lErjd zm8J}IOn=piQ$hYY7(v(!TEX7uIx{`jv~3<|j~noeKk^0I4eVsSS%>T==bQmb?fo#) zF1k;;74iA7wjC?aNI9sEtYo>yjOx-%-W8+%`vG2=_B{tw*Gxz6iP%y-yAh7dr4lh^TviPejkNbj@}L=CL%x6wTk}fQqzZyHIx1`w~+n0&A=Q@ z^$Mt@RW7|AhWJdZx0OD;9s#?fnDqa1M*8V2jSptS+u_-s9G^F~Xgpr^N}kzcQ~*xM zIZJ&nM0h@k2sLtL4uW$6$L9St&^Wc-)qguOISz0+4u~K8@)~^#e`C3L6vMV@2-Prz3FSD8@NBJ0E|XoW}hiQ_|GNlcHgGJ z>aH9*!6$eH`R9tn?Y|iwU*WQIU%rh)GP^aL|gmpfOGikb&bibCIfvW#N z%JMdz2fiGMa69-7-raf0)N2U!?^8;0wDHjrI7jv%g`Gj|wV!(**JeY6vK}<0LYL9L zbYSxNv4W-=aCHxJLG%Q|>Bu*#d=k5d{gf$@kvNu!`}5b)cST7F<)BUBDw{GV!XJ60 z(?96`3`mN$yQW)##!JzUGbQ^izC)v>K>p5c@vVW!gpns@3R+Ppm+lWk@y-5B%%;6y zHuTn1_s)$%_A}iKxrLNRAyZ7ph~XS+*J?@t+sN)xD9o69?oT4($(KHf3Eif7IO=UB zxOlA_$M0ceBD}Y%6tFKku5_iM_%_6_xmpo63EUz-yN^Vn{M=f+cdth97(Q|~kwMUSZgI0D1mX7` zPz~-P&cpAq_kV_^qW!Tw?QYK=ejNwtV}~r07ZCmr1NAmRdx)Sg@+J5LQ9kC9?0Np1 z?kKRfOP?CscF)J5|K^_(k5sVtqWO_y%P1~6M>08=`i}xuhV`{fJ=DJuYvC7z>Waap z+6iWl3Z(z}z`+$itlpzGt5dSTJmS?(?(c8j*muH`G@of%td8uK9}M_k>Zc80?;h@C zinzUo#`S=~zEj7(lHngt|H$#zrMQ1e)>Gcw@C|~q-b_A59;m$zT#pVsIZlF0JehpA z%{p=Y?-3tYWoF00SBhfxz~@N+;z|C#PtUP>GM_YPJW|m5)HP|*;0#ue=wj1>E`tz+ z=V@Bk6jrYo2*D#1PJ6_0{4HPl6vsQS>(@ z73UD&9y9%MRG=opxTR0hl{4rZME{Zh`l#nPxR5F#dv67`_v@xvT;)jupuEh(0S(9S z_C`g$cYG!_49>C06>v=7N70pX@^Z(>|RoQ2M=xpeJ^ZY>RKJnclS4 zvjZSw;rXGdLx>Ni)Lw(>reruaCeZSlV zYun!6nj4zO?QdRuclh9qVeGq=^$KQy{3EH}_u6@O6q>zTNzryC;QIW-QXH0Qqd>WL z>CrP|)Ly-)zFrF4@BKkXmPIFm()+GKZ|l3L^U-#`KMW=jQs&m zIlT5F{Rqpp_vU)fz5cU6H37)OpTFmApz&)_nBG~MP6nHHoVLCkWq5o2c<+6ta2*GW zY_A*OdDLEhiokrc)mdQu>&Zc1>|X7b|DBx$YSm5;12@l)qH2mrp1fFNq=(h{SHDpp zcjW-;UklBW{eBo9D!$s2REfy`0Vj1o+UZfO&O@OggFEu)z-^87L)iX)wAJ-QM@nS> zUCHrUuZ>eGb&#)3+?u7M;j2oyKZLUAqeBgXZAq z-y|`tj{Uo|hj+?Rf6Z!n{iy!z1Ye||*p;(wfAz8ZXDsO&2MumSy5CsdvE_#Voogzu zpU#2u>l_y{SL*$x=(9O>HPyB z`IJ0E@-o8T=da#c`<)DfDak{U-6i<{O1_o(6aHlcly<68>uMlA?3?Dc2iW&Z4b4R> zvk%aCiJg2mu4p?BHk^|OgjJCL|J6VE)}uKKbeT_TJ7E34)vgCsLc=ZtQzq9-;7hxFX+;w~4%Eqy`SH;%`#y0J8b96l@Saa`^> zU9Ry1JNG|y!o;JUqYjrJv743-ebxu1N^6{Gk0QR+a27muqVIzfe{O=FZTA}oJLfJ= z>oVPt6r*@B#=Qc^Im|dtkV(SMBl;>Rbve=gz|3pxkr)37 z7<;TIaF1Lgj^})vLYRR>CcH#B@uFe@#gllI?$*UU^Ds6dFR8601-I`WBH=TP@p$DL z&s0&ZEpErGQf6HmyT5&fZ|VIfqZnL%%Ovlzf;S1?@7!N_mK&|tT91`kIFLKxx1gO3 zZWf8S{S`hvFIMawvY_JJ60Or;a65MrJXy6Wy5PykbityLQe5A&+3wHmYZA!cedeRL zDatqa3|$nc$#bx6ZoFxqavpA{;hkXLm49>4=09Bqm%nKKHCmM(Do*Q#@yEX!5dWaK z-Fo-Kt$+u!VADW8;fOop^FeCK32_G^bcn3<5iKgm@$5)yelIYGy)SboPu@Bg@!8h; z+?ka}voOad_Qd2B6gOVSEt#urj6qw?3)Fn&Xn&^MC%~n=WfT@Zoc4PG8gTrcJz@Np z7@7g^52~N045+_0^Y>3wVdoCN_ob%fe?@t9V3Zv92kiU3%4E;UM{h_tPPxpq*&Y%R zD#jeM`m17#%Vj*yiRobP7Ra|(+I|?T#Q(QxV>F~_qzkgvr3l|)Pr&7`l+4r1Z;$}( z?bn-UBRX(-m4cqK6XiSHki$32U1xf|DswAz&(__WvZ z-W+8;E`QdNcbg-U0Hmnac8OM@xa!TrFZ!Ko5Uvdyby18UoIMZZ7wN6AI^+Zi!@rM_ zKgU_^8A6-UBil=Ul_zp!hKRL*C>Z(<~H`jgp~x)rQ+YVzg63 zuYD9=9NwuZhpjue{x>{aTz_gu7tHK`sC`l?50}T^wb>|ZCj;eYA#b7s3UK*s6+z^Y z%?#|+UUOjiQUkZ_c&i^@eqk{TlROt9{tl!3$9ZooW3pTy+(=UVO!=h**B7X@6q>s^ z1NOxB={Pr`cs1!^`|GdPB&@J}8DK(0eCF0Uqm+K296C`Oe-Tead?vCLj^5ozge_G9 zpkL@KZl9}4mNMoo5tb=An9m(R`ysqyKQgIuW?&Ly>6iufuF)3$?Lp6s9ZC2v`)L z@pbjC~j2=JFxuTHFq2>>aDj$N_Lqg70>D3=x+n zo?P`^-Zc-KU&(XHn0^IY|0~eD`?Sb3_Pw6;qOv&w$+dF6lv@9sfwwQFxa^ujd`K#h zlbyiMv4Z^j8N#r0tt~s^k*ZUYWFj<4l^*xew!`HktDz54J;`u0oO)Me1sX@cQm^=^ zr}aVm{9`-xLP~L*-Y+Z%$)VGL=D&qpjhIPX{^8M?(ynLS*nR19 z_ae$?xVt)URGAZDS}Xrw-HR2te>9(P9F~If@QRul?{X;Gk2Vl24(SoBPmY7;s}Po=f4p@CBDu& z*s!!fONF2*?u$t4Vs+WK_$CGwiLrO);IiV=D)-q!9RIt(f42|x&BA)C zsyBc8$8dS6vx4w1Wdc0gqP9VIyd0NDI;g~JT*TUS@aCfAj%i%pNSHZik%`q|V{@G@ z4MqGsB>lLfEd=}SqeC;&KNj(4-HiIFYuq4wHQ+oLuYmfcWFIY9*BS#ZyZP&h;wXNK zH@DsVRWJ+PA~)VEy+Y@!c29iEK6SvZ z+k8YU9v;x!h`45s%Y8#f9co^aVO~z2i1JS~E@pnm`Ri}?!p@GW3mp__d>x|tn8DI8 z0}_Pn#fVZUKJzXahz5@mVbAL?iIEQqar+dP#F$>*$Lg8vwdcU$L1d%!4Pu*@r{l!|w(>_Zz4u}5f22)uf{lwCz z?+lm;eV#EcfRO__RV^r@tA1zrRHuy9nA8Y;xoO0#^YJBb1dbw zldvz~T8h(CWPjd-V?4C65C;FY>fWeFcrJvW39{hnhe~Zc_Tk4{@OFjJd~tk-?Pty2 zy(yt(gVu3fQ5P?k)6YW%<*P9yRy6P2CF;i0VCOXvd&rUG!W`Vr!KYq9hH{gTb4It( z;t|55EvM<$m|Y5YD~a=KVDG4G^{*4H<;fz7X{_kCs`QT=^>9m0tN1u7`e-NY8ST_y z-gX~eCP(BZ;eJ1)U-x&P*G2eej%&ZzFq{I^0uLur`Vjt$Zs|%dy+@&*f#<&0%gFwI z{g@;Awyl8b*1EY+EsCFM7X!PFNq0j!Zr$fKktsNybI@R*GJPJr4fCgxA$8;OY^{T5 zJO_#JoT^M}m!>T)=YJbgapF4GPjzMcpSPmCt|Hj=_pS79So=KhWsd{mVc(pE)0HtI z$W|SeRP;snNhT&nmv1EWLWjO+7YnIk+&)RDI#lk-EO^FY*Y@ZsT5pDDoQIOtL|7QE zO-eEh!1c|DI$|mtBv{or6sAqpg3F~tFP$&n*#pVyX*Dfv)nLp2mUBAMZr6w)#gml8 zMjXWDPrcoe?AD0z^!MK|`J^{4S9>WNsr?0ekL)DP#a{n1+`dswSYTW!5ju{^zH?&; z#^wL98dwHRli+Wi3#FVN5gwbfhlFk|b-@cctAtByC~llLtgu}WA^}R7^S`a{BVKuY zFSb`3>4T@nUs|tIq5mCjkPu`(KMnGyh%*0CPT~K{5j2|NZ#WL`i?80Dk!i%`s+1nt zt&g&SdsD6?=R6u;udvTWA8L=mGFHAEqXIM@{W2oue5I;EpLfF^8hfPA=rdGgVAKQm zlzGy)^yTCBNw!erV#FNa*O)xkesTnti`|aN=3t+L9a2|XzrI21i>KBz`SBFkenNG^ zRKprtUmObkay=ib>sZsqAn7uP;ylv~`3(`PS1W$=yV||G1l)cF!}w_d<|$~iaqvnH z8}i#7y=$+|{MP`rh6--g2&3^hemsz3gRT!Y7+cVNpFnsF?t~u?{Wb&W_Vf&rZcl(M zo&-euP_iQsS=k z&BtQval_uT0SmFXJSe`;n9`aA8%ajzv_7N#hm7oh^o-bh)qi$mUGQ*E#r0KBB@s%( z=Yg=pUOq>4bYC(0K*L*Ks!4cqK`d&73-L$yNzhWURU!Pxd}`@(BjST;#M+Bzk)v?% zK=Cz}8I+IN3yZf}thRtVmw00*`bfART%DWLWf=+Zq^5V%+eH-L)QAVt{WkmAY0>^s_~*#&iflU6zfONmMvCm3KvCglf)6X|U(LDZsG!~R&^+K{u1GHG zuNx-$4r@zTy_C(HlzxS&xc$=^eGGgQ^RR=LJ^u*^#^tRhMy)m#WT@bFQ-iu0#U(+I zBK7zCEF2M|z54bo%9oaIMC_hEIRRH(ynfPmqI`+YHi*qKrx!|HRfQqXTgZAvhk>s701T*!M@q44o_rs_6S8_8opEbXo-H z^hAqdNE7*!_RaZ34$CR1uC^=fM=P2y)J`0_vHram2JWup^W#Q%N>AOPUcNF54hiiR zV~a<4xE!LT^+Jg70;`-E?RyljPB69X7r#M*gIA21T%%FHYoEQ>-n`Za%e0wWd;c}# z_%mcbMl@jO(^s<056GQDc>ME{vP^w@;2#D{(%IZvT%SYKXYcO&b0AHi{rAI|5nOI6 z!sFC;a1`or=!-;uMF0E#M~l!&+YS)5k)O8fHsTxKaHX7>0|7oyT~1K7L-AAULHE1d z-YF<9{-}Dz2*pos@z;wtCb4+ldvWz?!wtAK?pB)H1ELMF`^o?Pa)0*|t-CAS|B;+? zM_?!$+hC(I!jq6ScpqZ-!)}W(KxDLGq13&C&H{>B4&xP@9_4zxP;0m(UYNx5rfZ>gUz`7 zv|C%f05kRu=H;&RrNnj|r-^=<$M1*%sJ8!P^+F60ms^W+Mmn&K!IiL;uc=4S|IP`7 zDUE@8F!E3^MDhw6zjk(?eZII7p~fd`xjLnC+|E!^$L|AsNU$I-bg_LN*=c$v{qd6T zAna*VIiGwMt#=NWfPNq6F)*9usQZkIfZI9RF6teQ-3#^gtFn&C3B=_Z9n8mQAsJdH zzEf>rFM(VAE^o=#<=QX}X+__)&P|~9YDb@zTCV7W`U)=9zxSa2DmuZcYj$t~9yPsp zKG6^1luJ8tNkPv8?)e~p!wx$K-m)(|5xlq?`_3`ApQ@va1>yJc5lzw?9)sD<9G~O5 zXL0?J0_#rcuUOsM0jr;%vJoFvJ>{yihpM6Yd=%lvQ^bb`HO`km=x5+n1)!nx1>ukQ zaGyHv)c_PT4zn%!g!+qlI@~D*t4~xjb-wRy3i>`;)luN-Cm8}LF3Fr_zBz{bC(eP5 z`Dpn#bT88LuTw?yh4i(X3jL=e!7hD%-bt)3`&N4`6gee6VBe>mypk2n&x`O|@Qru% zXtltq!jiv*WR$Olue5g%vH9)~?H)qqPQ>S^(%JgW=UBcc^fScjGisMmkI1Q>u331D zEe;$}LGhWikU`6Xf4&!kb8tJ?s*kVO(oH~T{@-3#)lmN)G*=&fWv352rB1T0W9KAWc+4Zq zTY9nc4!1JDfDO?^+1tKKBb9*jC?ILO}l4*R&b9 za^w5p>?#^xx>AiJB@QF7?p)UU6Bctg&cIQtlCj`k`0k|()2DfaN5o^~O@Zqy*mv&` z<%|l#v-U_}F3!IT+y57Inwt#4aVpsNpSaaa0xu3R*^TaK#pSc7jlY~Nn}$@bHGl4t zw(<7w&8J6S^}(4x%TsZOQ9KcFxYqChPQa7lA{KEHsJ-gfV`pufa^XEehvms__fgGC zzvxoi^~174jf4{Ma{Pbam`~2uoSp)<2}vI6Ow+i$_g=mX@j4Nz+>mUy+LejR^I|Q< znXb%3<4XySAIH#sCeeWI-{&iHkYIG;NQ)2hkC7@1uSC);Or`v&pLqoNXY4pbh6E20 zYRNff%Ud8`&1$`s4aB~qKggY#rff5S|8KXck=rZoak%p1{;<|D)ZTX)r(VuK{{gsm z8+|qDMtnOaB6=kP48em`U114XsDBygMGWJxeUqiJBY)3h`zTxd)cYR%_*2~=WO{Ay zdsrXouM}2QUt=8u7X*)gQ7}gOJ)}csQ?y) z{r54}s|#A$IrM`iq`z)d{L`~%269uEaS%!o{@2em!auQ%!Vdn%q%tPNXT1vtohDd& zp!4CY77d}fI8I(c-(p!|A_(Q6v9xtPs}qJ6rWSS?R?j=zZvFM zefDv}@`A1Q!Wa5(NlX*)b_Flf>22}BtmTll;Ad;-yEj>Fa$9_uBCvbT?HY#cxg>_u z!ifL3bS&Z&N`}ELHNKYW8ic=^+x++Be^{NyVMCJz*Lbk?zgNu+EO%k=R;E>@I0pWg ziR0m{$dUiPM1pU$-*cPSIt?X zFYLw6Wme1PmKp1DJIXvNWy$xtK-t-uwvT`2arp!1x82@1hhf?zzt4%=$Ukwea&ppj zLqLD1>u&nC@7eFGCpGRnJORfoNCtk}@(t>LzpfNUdPDVNUJKjnl`$PX7kGSR62_Ye z^?&t4{ae4LeX&fV4i-O5IPIy5`uFEfY1+M4r=WFdoVn>b;-~W2SCUHDK8w$gTafAo znvbj!&FHi)_d@qtTM?DPdfY#w-^54!*a_gSI@yiJa}1X+d$1A)9ui^vD7SAxZv`%w zDx!sKy#U9JV?BEV)!8{sD>x7C?< zO~HYw94|gMcJ zdr#WC!)X$i&k>S%?e>s>fy?R0Zwx3e((ufCcF6*J&&YaA_B1uZb3KrF{97skJd>cK ze?Qr6_+lqwk*R+6{MKwI=~)Cw0qaKg6r4*;6}z-I)U+`;-IiLQ$T; z{aH!lA{K`uOU^BpV)x3oa2jyMzB+D70B`b%nw-xO{^o!rDm!i>G@oM*=Lkmegz$53 zX6>Z~P!#@$>{UOAw^zi+ZR!8+L*A`v&3t=V1h=olXt}7(vH*zV*=xe3sJ%MIe5bz1 z6TusY=K7R%v|s1+Hc#GznE>VjeI89zqIqPta&$ZztLJv~nW)u05yY#rKfQQ}tn=W# zz}@=hBdEXpS=9Ph&!)qDmz?dM7h5O;(W52!tYcgma z4)vr~Li#s8w`=nROoE&@pF{skK>116sP8nN$O7m{X2`fzfbb~m85}sTP6C4=W#`Z8 zqqs!vyScG9lLUPC^3iiFBmaxaU7UP6-w!mpcf?L_dw=0xGQH#Ce@zHDwua7wB-&z`8ErIgxs z%1kc-j5<596JO;2 zuJ4nDX@Omk7%E%rvxw$Lv(WwFLI>x-`W^Gfo@9jo-E{2-sdJN1oG*=XXlWY9lWpji z^V*95E_2+I9leU;NzW>CS?n|c=5qhmpt*z2*Z!Vwru%er4$#k?il2)__FV|dOL%{iQ^tFUjhcRc6uwKELaEuB3+I ziQ0Y(dvk0Z>Yw4DZ&-%%_V9S-bpJir`o&n(pnd}F?@8At?DN6me-7>X#)EDAuhQw! zzr{NaN9w=s|NaT_f9FNzOM_T@u?nw=Aj^HfBf&nN1x(NpXu zfJ2+o4}NE&c|o?$E#51K43B?*(G~K(6xVmDwXZI6CWFy)R93ZE{g^FYmHIPr)p1V& z(t%%oE{9Qk;MlP&l!W=G%xXixc$@PObGw66)W@%U#LGu(~Hgt%_9a)rdclI$aUF)QO;FMONlf37U_5 zOWFTi?c@`H0cx$bopukA{*JHq_k*!}Q3{k3=_lQyh#S>+P|II0k|>?WN1nVm+`%~uX=xC`Hloc{Lbj`y>({ML;Gq;-lsD+Bl>p?w}j))_r9`^5YCR6 z0bEZ63!lA$*01>YzS@64Xlria6okJ{P1&$A^)1l+HzGe84&j8+YTHgE9K@*|=ym46Rp-<23m#qdifP4|{ z0kxUVU(+gI*+hi6?_WzhB-*#0yNZ3{w(KUQ8H->LY)erllRl7DIl zaW3=sioPvx?N6y%_K593#QTOlUz5H5zOui3(uyP24B?(5+xV}?1Mfv({`u9g{v5(R E0Fb|Qp#T5? literal 0 HcmV?d00001 diff --git a/Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp b/Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp new file mode 100644 index 00000000000..324af538dce --- /dev/null +++ b/Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp @@ -0,0 +1,104 @@ +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include + +using Kernel = CGAL::Simple_cartesian; +using Point_3 = Kernel::Point_3; +using Vector_3 = Kernel::Vector_3; +using Point_2 = Kernel::Point_2; +using Vector_2 = Kernel::Vector_2; + +using Point_set_3 = CGAL::Point_set_3; +using Point_set_2 = CGAL::Point_set_3; +using Point_map = Point_set_2::Point_map; +using Normal_map = Point_set_2::Vector_map; + +namespace Shape_detection = CGAL::Shape_detection::Point_set; + +using Neighbor_query = Shape_detection::K_neighbor_query + ; +using Region_type = Shape_detection::Least_squares_circle_fit_region + ; +using Region_growing = CGAL::Shape_detection::Region_growing + ; + +int main (int argc, char** argv) +{ + std::ifstream ifile (argc > 1 ? argv[1] : "data/circle.ply"); + Point_set_3 points3; + ifile >> points3; + + std::cerr << points3.size() << " points read" << std::endl; + + // Input should have normals + assert (points.has_normal_map()); + + Point_set_2 points; + points.add_normal_map(); + for (Point_set_3::Index idx : points3) + { + const Point_3& p = points3.point(idx); + const Vector_3& n = points3.normal(idx); + points.insert (Point_2 (p.x(), p.y()), Vector_2 (n.x(), n.y())); + } + + // Default parameters for data/circles.ply + const std::size_t k = 12; + const double tolerance = 0.01; + const double max_angle = 10.; + const std::size_t min_region_size = 20; + + Neighbor_query neighbor_query(points, k, points.point_map()); + Region_type region_type(points, tolerance, max_angle, min_region_size, + points.point_map(), points.normal_map()); + Region_growing region_growing(points, neighbor_query, region_type); + + // Add maps to get colored output + Point_set_3::Property_map + red = points3.add_property_map("red", 0).first, + green = points3.add_property_map("green", 0).first, + blue = points3.add_property_map("blue", 0).first; + + CGAL::Random random; + + std::size_t nb_circles = 0; + CGAL::Real_timer timer; + timer.start(); + region_growing.detect + (boost::make_function_output_iterator + ([&](const std::vector& region) + { + // Assign a random color to each region + unsigned char r = (unsigned char)(random.get_int(64, 192)); + unsigned char g = (unsigned char)(random.get_int(64, 192)); + unsigned char b = (unsigned char)(random.get_int(64, 192)); + for (const std::size_t& idx : region) + { + red[idx] = r; + green[idx] = g; + blue[idx] = b; + } + ++ nb_circles; + })); + timer.stop(); + + std::cerr << nb_circles << " circles detected in " + << timer.time() << " seconds" << std::endl; + + // Save in colored_spheres.ply + std::ofstream out ("colored_circles.ply"); + CGAL::set_binary_mode (out); + out << points3; + + return EXIT_SUCCESS; +} diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h index 1c25058d84c..cdacadf1c5d 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h +++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h @@ -20,6 +20,7 @@ #include #include +#include #include #include diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_region.h new file mode 100644 index 00000000000..afd2f2c64c0 --- /dev/null +++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_region.h @@ -0,0 +1,340 @@ +// Copyright (c) 2020 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// +// Author(s) : Simon Giraudot +// + +#ifndef CGAL_SHAPE_DETECTION_REGION_GROWING_POINT_SET_LEAST_SQUARES_CIRCLE_FIT_REGION_H +#define CGAL_SHAPE_DETECTION_REGION_GROWING_POINT_SET_LEAST_SQUARES_CIRCLE_FIT_REGION_H + +#include + +#include + +#include +#include +#include +#include + +#include + +namespace CGAL { +namespace Shape_detection { +namespace Point_set { + +/*! + \ingroup PkgShapeDetectionRGOnPoints + + \brief Region type based on the quality of the least squares circle + fit applied to 2D points. + + This class fits a circle to chunks of points in a 2D point set and + controls the quality of this fit. If all quality conditions are + satisfied, the chunk is accepted as a valid region, otherwise + rejected. + + \tparam GeomTraits + must be a model of `Kernel`. + + \tparam InputRange + must be a model of `ConstRange` whose iterator type is `RandomAccessIterator`. + + \tparam PointMap + must be an `LvaluePropertyMap` whose key type is the value type of the input + range and value type is `Kernel::Point_3`. + + \tparam NormalMap + must be an `LvaluePropertyMap` whose key type is the value type of the input + range and value type is `Kernel::Vector_3`. + + \cgalModels `RegionType` +*/ +template +class Least_squares_circle_fit_region +{ + +public: + + /// \name Types + /// @{ + + /// \cond SKIP_IN_MANUAL + using Traits = GeomTraits; + using Input_range = InputRange; + using Point_map = PointMap; + using Normal_map = NormalMap; + /// \endcond + + /// Number type. + typedef typename GeomTraits::FT FT; + + /// \cond SKIP_IN_MANUAL + using Point_2 = typename Traits::Point_2; + using Vector_2 = typename Traits::Vector_2; + + using Squared_distance_2 = typename Traits::Compute_squared_distance_2; + + using Get_sqrt = internal::Get_sqrt; + using Sqrt = typename Get_sqrt::Sqrt; + +private: + + const Input_range& m_input_range; + + const FT m_distance_threshold; + const FT m_normal_threshold; + const std::size_t m_min_region_size; + + const Point_map m_point_map; + const Normal_map m_normal_map; + + const Squared_distance_2 m_squared_distance_2; + const Sqrt m_sqrt; + + Point_2 m_center; + FT m_radius; + +public: + + /// \endcond + + /// @} + + /// \name Initialization + /// @{ + + /*! + \brief initializes all internal data structures. + + \param input_range an instance of `InputRange` with 2D points and + corresponding 2D normal vectors + + \param distance_threshold the maximum distance from a point to a + circle. %Default is 1. + + \param angle_threshold the maximum accepted angle in degrees + between the normal of a point and the radius of a circle. %Default + is 25 degrees. + + \param min_region_size the minimum number of 2D points a region + must have. %Default is 3. + + \param point_map an instance of `PointMap` that maps an item from + `input_range` to `Kernel::Point_3` + + \param normal_map an instance of `NormalMap` that maps an item + from `input_range` to `Kernel::Vector_3` + + \param traits an instance of `GeomTraits`. + + \pre `input_range.size() > 0` + \pre `distance_threshold >= 0` + \pre `angle_threshold >= 0 && angle_threshold <= 90` + \pre `min_region_size > 0` + */ + Least_squares_circle_fit_region (const InputRange& input_range, + const FT distance_threshold = FT(1), + const FT angle_threshold = FT(25), + const std::size_t min_region_size = 3, + const PointMap point_map = PointMap(), + const NormalMap normal_map = NormalMap(), + const GeomTraits traits = GeomTraits()) + : m_input_range(input_range) + , m_distance_threshold(distance_threshold) + , m_normal_threshold(static_cast( + std::cos( + CGAL::to_double( + (angle_threshold * static_cast(CGAL_PI)) / FT(180))))), + m_min_region_size(min_region_size), + m_point_map(point_map), + m_normal_map(normal_map), + m_squared_distance_2(traits.compute_squared_distance_2_object()), + m_sqrt(Get_sqrt::sqrt_object(traits)) + { + CGAL_precondition(input_range.size() > 0); + CGAL_precondition(distance_threshold >= FT(0)); + CGAL_precondition(angle_threshold >= FT(0) && angle_threshold <= FT(90)); + CGAL_precondition(min_region_size > 0); + } + + /// @} + + /// \name Access + /// @{ + + /*! + \brief implements `RegionType::is_part_of_region()`. + + This function controls if a point with the index `query_index` is + within the `distance_threshold` from the corresponding circle and + if the angle between its normal and the circle radius is within + the `angle_threshold`. If both conditions are satisfied, it + returns `true`, otherwise `false`. + + \param query_index + index of the query point + + The first and third parameters are not used in this implementation. + + \return Boolean `true` or `false` + + \pre `query_index >= 0 && query_index < input_range.size()` + */ + bool is_part_of_region (const std::size_t, + const std::size_t query_index, + const std::vector& indices) const + { + CGAL_precondition(query_index < m_input_range.size()); + + // First, we need to integrate at least 6 points so that the + // computed circle means something + if (indices.size() < 6) + return true; + + const auto& key = *(m_input_range.begin() + query_index); + const Point_2& query_point = get(m_point_map, key); + Vector_2 normal = get(m_normal_map, key); + + FT distance_to_center = m_sqrt(m_squared_distance_2 (query_point, m_center)); + FT distance_to_circle = CGAL::abs (distance_to_center - m_radius); + + if (distance_to_circle > m_distance_threshold) + return false; + + normal = normal / m_sqrt (normal * normal); + + Vector_2 ray (m_center, query_point); + ray = ray / m_sqrt (ray * ray); + + if (CGAL::abs (normal * ray) < m_normal_threshold) + return false; + + return true; + } + + /*! + \brief implements `RegionType::is_valid_region()`. + + This function controls if the `region` contains at least + `min_region_size` points. + + \param region + indices of points included in the region + + \return Boolean `true` or `false` + */ + inline bool is_valid_region(const std::vector& region) const + { + return (region.size() >= m_min_region_size); + } + + /*! + \brief implements `RegionType::update()`. + + This function fits the least squares circle to all points from the `region`. + + \param region + indices of points included in the region + + \pre `region.size() > 0` + */ + void update(const std::vector& region) + { + CGAL_precondition(region.size() > 0); + + auto unary_function + = [&](const std::size_t& idx) -> const Point_2& + { + return get (m_point_map, *(m_input_range.begin() + idx)); + }; + + // Use bbox to compute diagonalization with smaller coordinates to + // avoid loss of precision when inverting large coordinates + CGAL::Bbox_2 bbox = CGAL::bbox_2 + (boost::make_transform_iterator + (region.begin(), unary_function), + boost::make_transform_iterator + (region.end(), unary_function)); + + using Diagonalize_traits = Eigen_diagonalize_traits; + using Covariance_matrix = typename Diagonalize_traits::Covariance_matrix; + using Vector = typename Diagonalize_traits::Vector; + using Matrix = typename Diagonalize_traits::Matrix; + + // Circle least squares fitting, + // Circle of center (a,b) and radius R + // Ri = sqrt((xi - a)^2 + (yi - b)^2) + // Minimize Sum(Ri^2 - R^2)^2 + // -> Minimize Sum(xi^2 + yi^2 − 2 a*xi − 2 b*yi + a^2 + b^2 − R^2)^2 + // let B=-2a ; C=-2b; D= a^2 + b^2 - R^2 + // let ri = x^2 + y^2 + // -> Minimize Sum(D + B*xi + C*yi + ri)^2 + // -> Minimize Sum(1 + B/D*xi + C/D*yi + ri/D)^2 + // -> system of linear equations + // -> diagonalize matrix + // -> center coordinates = -0.5 * eigenvector(1) / eigenvector(3) ; -0.5 * eigenvector(2) / eigenvector(3) + Covariance_matrix A + = { FT(0), FT(0), FT(0), FT(0), FT(0), + FT(0), FT(0), FT(0), FT(0), FT(0) }; + + A[0] = region.size(); + for (const std::size_t& idx : region) + { + const auto& key = *(m_input_range.begin() + idx); + const Point_2& p = get(m_point_map, key); + FT x = p.x() - bbox.xmin(); + FT y = p.y() - bbox.ymin(); + FT r = x*x + y*y; + A[1] += x; + A[2] += y; + A[3] += r; + A[4] += x * x; + A[5] += x * y; + A[6] += x * r; + A[7] += y * y; + A[8] += y * r; + A[9] += r * r; + } + + Vector eigenvalues = { FT(0), FT(0), FT(0), FT(0) }; + Matrix eigenvectors = { FT(0), FT(0), FT(0), FT(0), + FT(0), FT(0), FT(0), FT(0), + FT(0), FT(0), FT(0), FT(0), + FT(0), FT(0), FT(0), FT(0) }; + + Diagonalize_traits::diagonalize_selfadjoint_covariance_matrix + (A, eigenvalues, eigenvectors); + + m_center = Point_2 (bbox.xmin() - FT(0.5) * (eigenvectors[1] / eigenvectors[3]), + bbox.ymin() - FT(0.5) * (eigenvectors[2] / eigenvectors[3])); + + m_radius = FT(0); + for (const std::size_t& idx : region) + { + const auto& key = *(m_input_range.begin() + idx); + const Point_2& p = get(m_point_map, key); + m_radius += m_sqrt (m_squared_distance_2 (p, m_center)); + } + + m_radius /= region.size(); + } + + /// @} + +}; + +} // namespace Point_set +} // namespace Shape_detection +} // namespace CGAL + +#endif // CGAL_SHAPE_DETECTION_REGION_GROWING_POINT_SET_LEAST_SQUARES_CIRCLE_FIT_REGION_H diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_region.h index 94a318e46be..a142d34b789 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_region.h +++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_region.h @@ -264,6 +264,9 @@ public: { return get (m_point_map, *(m_input_range.begin() + idx)); }; + + // Use bbox to compute diagonalization with smaller coordinates to + // avoid loss of precision when inverting large coordinates CGAL::Bbox_3 bbox = CGAL::bbox_3 (boost::make_transform_iterator (region.begin(), unary_function), @@ -276,6 +279,7 @@ public: using Matrix = typename Diagonalize_traits::Matrix; // Sphere least squares fitting + // (see Least_square_circle_fit_region for details about computation) Covariance_matrix A = { FT(0), FT(0), FT(0), FT(0), FT(0), FT(0), FT(0), FT(0), FT(0), FT(0),