From f1ae4efb038cbf03aa807cc6597cd6257f354013 Mon Sep 17 00:00:00 2001 From: Yeicor Date: Thu, 30 Dec 2021 17:56:16 +0100 Subject: [PATCH 1/4] Testing face-based mesh importing --- .gitignore | 1 + examples/monkey_hat/Makefile | 2 + examples/monkey_hat/SHA1SUM | 2 + examples/monkey_hat/main.go | 54 +++++++++++ examples/monkey_hat/monkey.stl | Bin 0 -> 48434 bytes examples/pool/main.go | 19 ++++ go.mod | 2 + go.sum | 10 ++ obj/stl.go | 167 +++++++++++++++++++++++++++++++++ sdf/box.go | 42 +++++++++ 10 files changed, 299 insertions(+) create mode 100644 examples/monkey_hat/Makefile create mode 100644 examples/monkey_hat/SHA1SUM create mode 100644 examples/monkey_hat/main.go create mode 100644 examples/monkey_hat/monkey.stl create mode 100644 obj/stl.go diff --git a/.gitignore b/.gitignore index 2f3d9387b..3d1a2a1b4 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,7 @@ tools/parser.out tools/parsetab.py *.pyc examples/*/*.stl +!examples/monkey_hat/monkey.stl examples/*/*.dxf examples/*/*.png examples/*/*.svg diff --git a/examples/monkey_hat/Makefile b/examples/monkey_hat/Makefile new file mode 100644 index 000000000..d6053c1f9 --- /dev/null +++ b/examples/monkey_hat/Makefile @@ -0,0 +1,2 @@ +TOP = ../.. +include $(TOP)/mk/example.mk diff --git a/examples/monkey_hat/SHA1SUM b/examples/monkey_hat/SHA1SUM new file mode 100644 index 000000000..688b1056a --- /dev/null +++ b/examples/monkey_hat/SHA1SUM @@ -0,0 +1,2 @@ +6712d16e1047f5281e2968ee4543a9aaac46fb40 pool2.stl +e2d04358f90044d6b5f73cd17df3e9c55e70267c pool1.stl diff --git a/examples/monkey_hat/main.go b/examples/monkey_hat/main.go new file mode 100644 index 000000000..7cbeb1c96 --- /dev/null +++ b/examples/monkey_hat/main.go @@ -0,0 +1,54 @@ +//----------------------------------------------------------------------------- +/* + +Imported monkey model, with modifications + +*/ +//----------------------------------------------------------------------------- + +package main + +import ( + "github.com/deadsy/sdfx/obj" + "github.com/deadsy/sdfx/render" + "github.com/deadsy/sdfx/sdf" + "os" +) + +func main() { + + // MONKEY + // - Open the STL file + file, err := os.OpenFile("monkey.stl", os.O_RDONLY, 0400) + if err != nil { + file, err = os.OpenFile("examples/monkey_hat/monkey.stl", os.O_RDONLY, 0400) + if err != nil { + panic(err) + } + } + // - Create the SDF from the mesh + monkeyImported, err := obj.ImportSTL(file, 0.2, 5) + if err != nil { + panic(err) + } + + // HAT + hatHeight := 0.5 + hat, err := sdf.Cylinder3D(hatHeight, 0.6, 0) + if err != nil { + panic(err) + } + edge, err := sdf.Cylinder3D(hatHeight*0.4, 1, 0) + if err != nil { + panic(err) + } + edge = sdf.Transform3D(edge, sdf.Translate3d(sdf.V3{Z: -hatHeight / 2})) + fullHat := sdf.Union3D(hat, edge) + + // Union + fullHat = sdf.Transform3D(fullHat, sdf.Translate3d(sdf.V3{Y: 0.15, Z: 1})) + monkeyHat := sdf.Union3D(monkeyImported, fullHat) + + render.ToSTL(monkeyHat, 256, "monkey-out.stl", &render.MarchingCubesOctree{}) + //render.ToSTL(monkeyHat, 16, "monkey-out.stl", dc.NewDualContouringDefault()) +} diff --git a/examples/monkey_hat/monkey.stl b/examples/monkey_hat/monkey.stl new file mode 100644 index 0000000000000000000000000000000000000000..b18dad2062041560578dec6b9b6ded96bad9bd41 GIT binary patch literal 48434 zcmbWgdz_YY{>MK;qGTjVsF5&|3^|Nyu4^`_oZ6AL9SI|$O=>nt&Lc@nYzR3I!eG!6 zTFrbuW4{hLEFsn*cF;P8K}&YS@A+!(_kF$V*Z%(a-H*q|^}0XT`}KN#&i6Ibc;AOj zubw<*>Z#+3PMb3M%%VXPPn|UG)G7P)>s!*dq~QP0>-kL!3hq61uk^n!mBu9({x%d8 zv?K)ugJMC!bgl(yK|u+)sjQ%&5ZotM4qby@opqq~Z8{2z#I%#(JW#?NPWjK^$s&jz}+%tOAn z+_Mq6?$N>N)fab&g*VlSfv)bz!;g>6lUC{zg_nkYJGG1LA0D1&>*q7>*nB9;_v!D8 zE=r#Mev{a(&#_{J)_GpX0Bb(PP%(N$QWMwF3g>MgBa-dp=c#~-6HDj zOf3FgNqXpre(~WUDN!vHJDt6*A47MxixC_HUF(zg{&0)mryp+dkEG8xkHGL+t9d^p z*Fx{T!}h~sVM+V6(+#(ZfrriDYOpa@g0&XviRpRb$mFf3HcPpt>(A{7} z#{e5+d=5)#+<>8R(CeF}i@QHg?Z!Nch#F53JI6>fhP78w-p_l`^N{N~8rI%o7&cZd zoLkl>DY>zNVdJByIia;rVta39F3HWH7~Ivkx8rJaE?w{E=A6itp_Q>NDu>1P-nmmW z(}+?JG2pBzhW?nkX36Yt@@VFbr;S{r(3acZXmku%rse0^_oM4+7dKMk;J%;ntmfN0KIfrv^L|8^Wel%7 z9zR>*-s6Y~7dKMkb#E^pxW0$!*5l#*h%U<*Uf2A@@E#9;@>1E(@17SYb=psNC^l9$ zK26klKNyV}V?*8i=@T-g`_H^3+2z3M`03SqZQx;ZzOGf=Gd7Ii7--oYoAV9(anMIE zC*Ag(5Gkwx-;dSkZn!_hh7mG`*IkXDG5lElYfby)%kRGmcdkAuxA%H-okY>AGuI0! zrD%8LI-T-r^FJCdS=Bba+WdRo6{EW`hS7!7E(<2p`ku=#lp z7%^jbj|M!9T?e*SJ-=i~I^^o2xZ7iwr`JcGODr)N1mW9cWPRqyPV9(~33ao#`Q zBEoz0PYVjV5v9I$2Rla07;2eCc|Y&*!B_p#Zd-PZ$Dgwlk1=RWv`KgWp)mgP-R@~Y!CGo`w;Q;GVo~1vDB*r>G97e75P1>7I?-#Z#wsa|KbjMy?>cGJ-K(5d$X7?eE2dHAvF@$)Ikk70ED<+m+j zpka(!aDH^bHf}IN#;Ap5QQpscJoR|T#)~ejh;Izj=TS3SwtjS;Yur|I)!3Eh?$^BY zURljM$46pS6C49<4DWq*_r;B`_uM76EFA{J?+LG~{h|}R39q83ZY!HOrXnsqL+6#> zYb)@umU!<2zC5Yy=FiTL7yM}-wE`{Msb04i^?RT=ACIBcD;tV0E{M+_rB6<^JI?@H zOT71jZto^NyPO_3KlEH?tbdxFI9|7idU=%NvE^U;gw2ke9Gi^uZiC@JGrX?mL7#^;JWF1<$vt7u{a7*Qhwz&vzmF^COhdDPl4pReCEokJy4}L& zyNr*Qy!$*lUQ0VwJuXMrd&5XFhL6&ou$Fl5&F1z^lel;6aLDyU`13=$gYmePtKaAI z-YyG6`yUstZ1`CWdTb|-bd4czJ4^0*adO#F&tDWLuDF9}ud66AG=}Xc+W*2Z8Xdz$ zd5?DZ8DqrTZ!ekL?%a5DpEIddd&|)(2-~>r_rT=FJ6IeyW*H?`hlHo7_)G;_^( z>vK@25x?~a`VGkUyS^ZX)1A$GV|b4RC)6b;-tfoRa?Ex_YfjaYm15^~b@J3uud7x+ zjF2(pA$Q%agk|g+urXeoaZl24;vXYdqJ4Y4uJobyv}y?3vJLPOL1i{sS; zj-lOpX&2Seidx=l#i&mU!K$AJ*6W=~7snn0W^c%7%ovI?_MF~VFPS>^;`s5<%ZTv0 ziqa_brTBiRKIb9E2D-*jlz9j<-stG%$4c)#LsbnXWB9S!cR*#j zvdck@dtcK|QAw-&c*dluto8RWb2b<&7%^is;h|4ZMQNX_?P4#a-bVM@{m)9D*tKW! z-m8B~TF)!p=}J^P-s`hn4DSsiW{h^+wWwt|U9r92ly|QxJLRA6C2MA!mhx$nYr7(r zfz1!jI$TKo6IU*4JhuCDNuS#$h=H#9v5w(9w8m7gHe_=lte09cHpk@qSpiaAr7v4?Q#Z zahIPl{5%crab2=-ufcI}4?W5Fc@|?hb+#vreQ1yC!ox5QZDr)fWZj%QMlJPX*U+ke zuHH3%HF%pibgL~>sb)PdFKhi2)eQQ1f{|v7a54M zyq5YTb`c)UVD&>Y2%|A$XiO~1>o!wwj31Ymr+=HVLp=S7v_qs!V{F+H|fGZ}rg*WAm&exy~eCy(<_-Ldhr!RN-(uBVBCR-OSihCFOU z9-X^+Xn)+Wc>WtZp`n*)`F?nB7)i$P{qSt>J^Rv4V&OJp;^T|1rQLdIZ)-p6{7_4D zLg;1B3@}p1Aks(W*^V)QpSrB9e|{9>p-;8x8~72OW{hI)TGWQ;m3cQl+ddum`N8r0 zzpNnI@AC~k>_ok`<^FJ9td2`M-bl~=KKFj5IqzDoi~n#>(x%sxSl0Ddo|iQT(pt%# z(xkP8it^qo_Lv*q?|63Hss3`JH6M8fk?NBel{qY>%SH`M`(D3w{Iz6hQcy6H($S*8 z_5@H1?nHfa;9q76t73g&bC~(j2x@0ilV^a9;rlV6%W>(1F7Je8%kF}q=O{61i1z*P z-j_apYV32@3*nl5j}n8(JOga~@Z)xN<8f*Ki1xAuFUxIK~F;HHN%}M^lzg+dY1L->C5I>jTpDJv3JO`J;lO zJ~>93F`CfIM}gH6ji~on|G7lEHoC%k($_BbelSuS*VfOdd=wseHrTxRTcU)&5Ne$E z*∨%U0Zm9=3Y)*GK<@fv)*UsTf6lPL)>uW#G{(&x-pUa7%7~*qw6vVdLZdU^HZm z^~WUB)fk#z(G^=Ioi08z4)5_O9$srh59>LM;23U9yr1`Y;-PQCptnZHwZAzPk9JA+ z6K^?KyR?XF2}QlG;tC_l7&>(Qt*(*ZxG1KCw?LOdaFW~NAqHa|b z#iDdR3eV>XBgq(Es~}`-ox8?3@U4BrJNsS`AM4b)VSNv^Mo~>%bv~-K%en6M{>)_f z%V)>UuV2oLd0l;G7{hy9_v-uM&g}=q%l`Be5juBu25Mf)!8Qu=&>Ak}x^m~e(_w$z zK0Y<{%K(FtXMl|{9xT6Rln!6ITe$hie}z-~olK4QQ2WKMA!F_MhYgqB5lKkw1} z!1Ku_T@Q?-KJS8_bPtd&bF&((y|N4s?*}7fjMaGLqrk@Y-jCn7XFBNF9pimJeL=)p zbo=Dm-Z^dkd2<+b-Va917@djkOHtHX6pZb?``2%owtMq;vCpV48d<04?g2w{(FLs2 zXCB!5!H5|{y|<_$u(7@O-ki{+WLroZR-8F7est4i|$O&7c_d6zy8oy-jh9kTL2*7UliC$NP2V$u?~! z#^=X9giakrD^YVihaz3+dp{T+vvZ<-mBs7oE84 zJ?5fdjB2p_7Q=q@q??l~|9pDnN^LzgHV58A^*T$IwK77+@Vb6{gpKj*vZk=~hF3!U zx%#(M@3B0myOgWuBtPeaFDYxOniBtg{hjD|-Q~QOt{C29zdP&7y6@0EPDnNOUaMW& zx`XBG*YRn87@mB0U}L!Upz&g$YwqoMuRSk@_j_|`W7u}s*!bPrYl!yy&U& z4@QzPysqE#dd@aR*Vq3P8g`l%Z#Z-!XNmXdkki#V=nS?~X~>GtlffsP5|ifNvd_J4 zhn%h$-s9$j{v0m+;OJO#Xj^nvbJtE$@3D&d15nWOFpMN)c-_?@W3K`mLenbY|&!M?XUeYj^3hF8=57jPX<83xv;+O|K|hiT^ziwo}CTnu#0eLC0}-tUI0Uz5|X-8WV|xCM-I zw6sc!xO$HY-fLBsaD@>vhSx0*8M^{(jGI2%D&6$y!7;t%A!;>SYOne;jG|V}=USb^ zi@Cyx8Dm(;qWqq(&eeRVZ`s@;4xost5t+7y46<0q#-Va917=BDjavoxP?{m6+5c)j4Lo(>HQ&Q$v zQ;f$7qBOrdgSBq;(kI?BFpMN)G}C(@Wo*wl{g>n7npI1ZJsyehwKmoBwP3ZqJGghc zz8{PRUcOGH%~-hjih&TZlnzo&mO2>9o*!QcTQce`}N8F{^$0 zp96YR>pe6l%Zbu!2-q=V#?XqkDDUSzw%MmK+5h;($zP5-p6H%vZN;^SRoF}mTJUHB zdp{TA0$3bPkZ$;t-^%AA0tBhLa&A3;oui~-yF8v zacX+<)ko_mFJ5^D5qSpK{K{V1x8bsyW7EY46mFoU+0sZ%$<=E%&jVk-|D|E8VH48F ze|rZtdYOAM?~PFiHosXv?VMaSdvwxmlXKvgCars7`Ubk*4@Srs@-TOePe+(yd%xzb zH%W$9O;h1c9Ug7Ds_kA^^)I6^fqEDrV|cAPB8)9;jAJgZPB(w-mSp<2yQ0%TX@|}@ zir1Y-{WGw%%V5Ne;k7h2#-0Z@#)+4mnSNbzZ*pPzKy(&Sw08EVc-yRCjBThLZ{jmE}==-rMD#L{ zHj~KZlxpxYu8S!B8XChHvwEhdT)ht?=i6J#U1RwEb-3<~wEdr}B<^y-tgK8*?M^J^~kcP$HH=awVW%KN_;gVLKj8!{GBzZLaexgP%E z{A7=Je@I5XJ0fy|CfOb7F5>+&MCoOIA9c7dY;otuNzv8%3DRpR${6!uX-(8}eRRxC z$&*DLlF5thcdr&zNrZONQi@)hi&YdDA!BIVeU!018HxsUG? z&id1Lyz?Gayf=pL--~xoNI!b3FkChx65*de1~OxFSrz4!(cmhs)k#4CjF>U>9BNVi zS!6Wz#%SF7#$@Y}Q`5RV`=hZ6jRVkY=Bk&iC1Lp#QcgYZtuTG3lhgaYJDpZ(hBn6V z>+SBcHp$qtPY#X$I5i&VT`|D_;Ygdu^*N#xN{61eJN{cnflNlTn}Cb*7KuQqUvFU zjL|=4QGLP2?g=)=C%avm%sja?E^Dv9j`mtd@xB0E&2uk`R*7a8M#va?Znr3}+lzW* z-1KDE^x~;Ig!^0c`$VrL5B=^!W2pLJ^r4(9jF>Uxk&gnKt~`t}@|1(pE^`JYM}6G~ zzSok+Vsv$~seUC`V+kW>40+_Ez@{q?V_f^s1JgZYN!oOGdogp47BnLu=g{Mwh#7AueoTybkQHqgyGkcZ@c$~5i^Ew@2ZgPNAJB2 zzX_i3MMd(}V~5bz&+vJMaw6D|oA-keGlm~`zgN8X(YN$Zt4lTsHHRJz#iUod$Le3U@-S6ER25?jA}$uZ;J1Y$Jx#^$g#q348YmW4Ag#o$==1 z*>BC1ozaq4J*6wS1SQ`O7)i#c$HStUISGudPfY9i9nW2#+~5AN`01c+(efS}_R0DI zBV-KU5AWf7zVlhThX?yikN?=djCSkQjkat4y`SpM14gnw!|rV8VM^BxI5FLP)Bhx& z{rw5XN&C5pUKZow&t2_#wOHO0cbFWv`01)J{}cU9(GXgcXMiIbI9E z@~{)t`@x7A!=Jlj(9N^K*22FYeSUgpQ+aaX6KnX1URsOM*I7G)QUR``XbkO4{QiTB z;(p(KIqhZ-N^{|8uKw`B^!&4? zg)`gkhKJYHXxbcj4;V3Hc&$o2jN!e%+T!x$@kyiOYn7i6;rmg6M*+Si)XPsT98ZA} zGKO^RP6gN)zE7jB?2?}L=QKU1ME`GwJ||WrM6!;%gS${RT)VrsZyUSpdw9C=k9A_u ziaY~sj5XZr0K18QB`^AIOHy#2`#$+~JpYp8!%;i8Nq3#CpVc(qc}7k_ zqu@QhnYn-Z<$G5p50&mk#6tA7$GgJxe)Flfop#@xfpM?r{+w+2r(MJ#GS2{;pZ8dB z%`eHm$E0bOhdw02kIx!(ZA85{jF2(Bhu8JvcJ_Z7k`E?~jc=Fe?_S5zhib0NC^NwG zxGtw$3a;h4nsU-E_l74I{5E#%^AqpJ(U&{}Yz&mO>CoFA`0eG}%^dpwlW72__jMnnBtFn?q+=8|(`r$G;Le%ys_{tj64uC?p^?vK}p zui8zGHAmhehSvfc zg)6?6&W2;vfQMe*8%B~b)OL&V*ST%J!07jI})KiDE);LuRDW!+rj5O-5ly_`^RurTQMkk z26!nhd2h0Bc{=8& z^OJqQ-QZadGFU6PfT`co*Ym6O8;l^5VQ)o76yKo<#hepePa6Q zX}>ooC2t+1IUkL#J@5K{jG?}qJ}p2EM$8!f(Y3aZ1{-?}*cd~e+pX;2E34xJZ61d2 zpFby}TSc8l)t~t=8Z(AJ!$+W7L3EuM;A*bRD81KSnSB5Gf>6K3cr^WfuHgL~B5D{> zKW;EW#;8EcqWt*yar^D)eRVsNs5D*N2j#HX8aHG7?f$LP|D4z(a-~*38&smD zQT2YKs28Jza_3E_#iy6-7*4)>S!n&)AkWC@%ER~X>B~E%|LQgX-G{DGt$ldRG@Aq*$Hofbj(&W$~?b7vsR`Y&-O~Z&8!+VUuubgqxF7oSm z#X%F(-W~5RTlr!Wt@8WMkJ}jhjs4%&V`T4kG*u`n(<@3rY z>;vynimv8ZCyubKP#7^|cr9ULc<)zdl%zM;UKj3N_b=6c4m9WLbv0!wcp+DfI*gbx zH0~DV=VL1M#u#2WEM587qQ=vA=}xU*qu#@>G5vq9U-AFHw*KF;e*JjA%AB|Qb=GtJ z(!q~?)5vuLkHKkc3$@xc|7vZH5i`c%oR*@zpZC~u&Cg|READUDa?OS4X#Lpk^d5dc zzRUj;{_dO~%D9SwN1g#zOXR1sOSl$%{H#Bw$1HrSv23$DsZo2y(BBB_e6%R3 z<6R%`61vp>MGUm^45zDeUTi6SwoQ-p!0neee)NVuZM=uh#C1`utI>vA88Kt%+_fn0 z=RNFT>+*C@^Bim~yRhB~Ri2pPj` z-GG9z=Yx%L&d*20Yx~@uJTj#%wYQ_S5}g`~*S(YawbV_gz(_NO*SZ}AW8VojM&Y5q zjsJbIUE2JF_})X$qU&11>XV)PeYt7>Q+okFJc`0KF}GZb<|A#NU&C- zcJ2M&?H8Wi;eyz1pYH4yzYla?8KW9p52FUg9-Cf}oZI`H_`B&3^R5~#yQ7n05m7p$ zyx+}_yb=EI-SKhQaeIm3w5EZLQ3GDrKiC&x1AX0{PHaCH0QQEy{`9zF)(BJ_WJ$k{my89I(fZmVtk?J z_KoY0l|M@~ADVNGRXF;cF#f*NgXSbclSm~_IGn}7aW378kt3WTmMk~{dJvfJ(P;Gq~!+VTGMSfc0 z{_}D2i9dy3h8M-Q@9l)nSjudobe{XYQjzmeefN9kC*AMu8t;8I5hFSVSRVRZ@@u+% z<=$z(S*4M}-C3~gH1K=Hd&5w4j^VXt;bCm={cU}7GNZIAzIE-JeG0=U6)@c`c4~?bm!DoAPPW!GsIv#OVI~Z5f53O*mYOm|h_Uhlw55x93 zAzt#!t74#)XMoqh)K3uJ`=1>O(tgvAiJNbC4R`%#hWGOxFk;5=9{#h$pY2Uo&rd$x zYEb<3iq7=G??J!syf=)HF}$wdcYZ&9G(42`8~JqErQt$)Dvvt)qu*H_k;4Uyi0apB z6sU*Mm@$rkCC0j_92Q$1#yEAemZase?ZWdroJ8${H2Ync*R7*oKRZ-0(lA2C@LCH{ zSb?9S>cGZ$cKD=p#__Yls9k2zUhnY;ngx{ku>9}(X6|)n?ET#};r)M}FNV{71Z)iN z@vp}!)3Gm}RW`YL3-r8JE$=%}TJX^4yjFf2uKU&=8?U_clmY!a6!Omhp0YM(4DUB$ zPEBZd@c87nck5qYz1F3?*Gg=r2-k95O)?@Z3QeWDmKbJz2gF}&a3-up4! z`_p!1zr1=946ikY_u5@5G|W* zulF9jWc~HyJ$#>5ZP_8+?yynmYww=H>E-9#d-(M>zIT`K#laV*e=qLMN#?cuTJmcQ zMv^i7-1@eA?^V5DO$N4^kTx{EgNL7IzplI=jF2(Bho5IZAFn+)Jl;2K){=&8cVKi! zP%0R6?J7T3L#Wr>&!*gXFaLsd$Zr}?d~zQ#qGNz!utQ7 z>Gydr)U>+_m;o3eWB73&g08is7uXnm9+=bko3rmL8~mq>sO^eYv9$2>9_mK{^@Y^K zD9ae$qbmx=mY*@w_xDc!_Q^r3Ok>r|6ffzw+}3N6o8Gb{oa{+okTEV#^^w;-SK{f)SJ>5Yu`y8 zzWco};rk=0^}2eZDODZ1!rj4oriKwR#v)o~cZ7}M{T}aqUvlMk!{euSD!Z5;AL$N5 zw;fmS2P0$*uPbZ}KW?x5G_&y^laGmO59@@U*VQxVQPgP`D=3NjL7!lhWsC}R?amyq zF}&X{edmOW&-*bv_5R45hxOL$7U8i6_cb3dl8jM=uHEr|@-QX7bs@TMT|i5%(wx_# zSDeG%gKu_7-8Va+Weo4(y|=sY>(FcJ331ws3DJIhOyY85-7~U_J{oB5EfxBvz7&bmyYku7L zW{1>$vqQ)jzE8D{*+}!P3&DNs0(zRM&$uq4bO&#iuJ6Y;J0$L#9q=%Q?~~1kTFAF9 zG`eqHphhcvF;{uy*W?O$qE8+5iEnm@?wcJ##%KXs3s-}U;d1uV zR*idp*1*HjZQ|?#^75QqWjhbBJIQ(50;0<$NQ1jFm>PTfR-`l9a=4hQc@%+Rq#47sbjLiz@aL8LP7B&>W2GpYn|3*@p0{w`yRo9i#zc2) zti+fTT3ah6_uU*~P@@r5R5L|uOi?!HYOR}d?yd~24BfFg7sLC#`G=9|n}1lIbm*hM zq4&BiMCMl_*j6Hpm@!tvvO6tcV|YKZ%9pvU@@Vy>b(1K1_2N1uhqY@jppdyfxXktF z8lx24k%;LVFvu#OxUBLyj|oJr8%eJ6E2ZApWUdb`bA8T3j5=^zu4`$oy>x?>p;{9O6jNx@_&^3nU+?2k!c6hpfpDwY>K3B2&YI7c%u}bhk=&uBO zKNvA%c-=}OjZq7hhm?-_X;SjW@!Q88cOObkMbdhnmvg5Xtlk!aUk2+8F6UDkM#va1 zbJy-j&lnY8eKJUitn!J=DxdP#RNN^-t2?*|tmiiGLFW44GS?Gf4DaE+pZ{{-bm*r$ z#6xd-3&)=58e=3_H$Bm${yZt`sr6hxaC{eCo2wqpR^*1EU12$mJXOk-0v&%=JVVLq9v4pPpi@ zg=CdaT~>K%bGKXCI_ElrJAl0(nd_s=Tpu!q)|B~4#aejB)w`sWck{1 z4#J}rJPO`+bXnzDr{(CXKg+@4F&s(_qw);Mjxjwqg^|aR*#bAAwXeD}&zP~*wxxCNu z@r8H3A)v}&j*O$4>^=KJm z1-O(}HKXGF$SU9HvdX8~s#av44K_BJ>&sl`dNgexbOZOMG!gCn$SOauILj)JhwVu2bZ~?Xk%!PE2K&!X<==GnP= zvB_K?T;}?mt{6?=1LVPP^t>NgXaarZjuoFOR zJ4*`i8H~rRT*+J?UFLdvYz%n}Axita0Bnqb!)GNw4LmLWx6@ZMgJR_kTANxo z`*R&j(Mz8)<8rtc*C`Yj12cxA?2c#X3^t`F_Mesv>V0%PX#WoI^;3jCG4+|}{WSkQ zsnfGBjF2%j2X@Enw%}n(`Yt_}Ro;12=X58aI~%OCL~DS|^}%JXcUsM0eb(w^n~k?I zdYv^lxnl0w@wNSD@w}_Cm-|X0YVp%Z3#)%(YwTczjG^}0om#LlDq)$@AFjGDId^CB zk=@dX+UA_rx|SqaOKQN<*FJ<1GKM}U?2d1{?~~?f<5}f($6Dq4MCST|O&iTzzmR@t zpU=bB_sMw>;kCSnutuaWh0OH>T;_T=K7QOhL)c}NcVp$}U_%D04mZ6v922!fJ4r8V z6>Bh;}7H^DnlT-Va917+zO>GR88nDUnsa%w?6&>8cgpL(e|OAai}A%Uqw+ z6+@x#v&t7aKhGvdeByG%=iIfH{D_L7J}sc0JoLfkp{G@9i5Pkb*X6J= z$Pu5o9Pzo9Xr^kx!(r6q7{XE_4}EZX=<$%27NDRQ0X%T%Kc$ zgssn~a{Ma5nps6=jB-5ej@ZWVN)vzGGCAdrPO;~&!}0JV>UEoP9z}S-2pK~wQ69SE z{fekpp9-iaM||pX#5=8Kuv#^qD8CZPLmyopdgsxc(_OcAVcE)mvElV({?@mLx|zqL z=e4{?Eq*pPFha(tMc3}gBi}0LL5}z`mm?l6i}K@BPQ8sfdFUHm9(t$i=cAl@V`z`x z^1)R0`eiu#C}CR}z8^3`#_;{{46pQmgD+_ubLO_OaT|SN>6PDiv|U(hsUT@>mG&Zx zGRNSq+N(Q$r+TG-&lr$?e1G>ix@;~w{sfS&ol(VLeMU`35k|}y#njjxts!H0KXSw; zE=PQ>?b>UlNh`Y+>^;asA6*`L+MnO&&cl0?BR+9C;_*;RRub($sl6Y0=!45ckB2c9 zI=@mB)IYV59Pz2k5uddFhE8|9uJx5i^D#_X4zS zjp>vt#Jij;IpULz=ZM$p(-V!JtF<5XT&me64}I8p9(wI_eIkq11UAN_75&1$ubvjq z?{^H2`qbAbj-oV!^9-ds-cS7yL;Dm) zS;mlu-SN6L)SJ?#t9pdZKbjmjU3EvI(aD~?yfstg4aDP;hF7}k;tugzzN@g4B_kB%Spc!CJ8Ta;$M;&H-Da+DQ?wFa=IZ^(4If;N_>eKWIln@%wQ$wv z4nYl`x(_2{jQWt>sRJ8h1=y56zOQq#dCv*)zVfHx zd)*q|OSgsV9E$p+zQG6?V@}BK)PRlQ{l5Rh4+A^?cY6HO?i1mA-D2KLcRANeitSn$ z4H=^{WOs_e#_)dRZ%}yhlB1rt}Oy8@T-K&O?4Wsh4rpcnL4%N{;x%<%kcN z*4tpUO1)hUHU|0IgUjE}-Fy_-qS~dcyIH>Eh;MW`;?XjOU$5=((*76Y`?nqpn;i6q zc<=?vu;KT!-xJ;gMv^hSmiO>`ko@hzx9!>bjV??FLJ~uE=N2O#_(rX z2kONZ=f;cuAO$?1IECN;kCTS>X0dszrE4rZ+9N6zz5J0 zJ(KD*(i6u*FgfA}x*YLNOXs0o30?~L!`(aSTyj`2S-zI%^yJha}5V1D22i1e<8wz2pTeZPTbNV*!E(j1UxVRuTv z#_)bq&p15YzolDz;*%d4MXy_qp7pi@EG6BC5i>>wx^~BpkN10`v1c6n^IwxWC-hDA z{Mb5!E6~!DhS+Ptn#G;LtH3bQjL}T*?M@BY7+z_Q=Z=fv!KahQK2?TTudCV89F*d% zvFZ-)4Th0s42_lDDFPeA`xX6M6S8xf~dypeOaXI4AHHP=_-sEpj zUH*1-bux5AcMZ5BSZ9j&BS(C6IpXm!#scS854IMPzddpJ+tDjRSNfWtwM2df_I~7u z4=zVM(Z*1;`KH4pquX#I@~jF>UzQEPWJ z$HrI$HlU@ZeJ*WCn(xqj)KWBxYX71f*4)&9O$kQK7&S!N9ckMf)aJ01$lsp2{Ot+tq;@W) z65PUdJf#|})htEuJT=zC3UDD&n@t$!i35WE$*t5kUP0S^O{BR;qs z@zfYYQTmjdTPBn-$lsp0{OwM6K6o}oqpa5=@H1d?#0Qrno|=ud#50~Mnw4B$@>AK| zE@#JeQh6=!p}8qWO{)|}$QYH>*q!;{jkW5#@i!$8O>3VlXt1pEnz2PfMBc7U>l;xCa$})=9aTUcFRd4T`zCUtp6-BP&3K;DyP7R8AH#ZcE>-X458kX$lsp0{OxETfS%TBJ!NAn zh8*$1<%ka%L*p|Vb|FPCTVv#JFLU|Z(N&KX`59QNa8iy@3%&qMj`-km#1m}{>CWT* zX&W$VrvF&>@8gChQ%W|L@!dS+_YCDk*f1J1#%E~S9bse0&sufjGquTQdksqO>DMXe z-hr!H;&oS{Z(|4}WQru1FuA!T1Z)-67E_mSx6S#&O1^6fd6*LU+f$dn9gQxOI`rhXlB?#(dypeOx*YK_ zV`$%LXU(Q~?@rI&m6U$`oAlOAwXXd7`HV;u#DPVX0v_)=ryb%kvQ%fn_6M$8!UushN< zhWESa(314E`*w-Dch`9HHM<}>|>{rF|I7_ zA70t)ytHf07wj?ryO&{jbcdnyy&l{ZTmpuXWQ=+|?9MQ-F`C)QJHyZjEbVf6sGpZ6 z2bWEO?>*LVH{WxxwH8K_G5nbLvGN||R!?1Sb$ps=O);9I!CEoZV1IUzzdgGA?X=Yx z+IR9?cZ!15N8z(yIX9g$x-@jz_7}#=>yE*roVo=>X`dIs(3uD$W(?`t9sQoc7^A_~ z!dYFWrM<42l-#qV3tC=x1bP+J)#W_w`F`7F2{F&(e z$lo4a{&sYY;m<4YM{f1h~3f+ zhozTJ=@xJkYs>coM$8z#z0E`$qZ?|bbkcseC7sSZGS1qrEwx@(<0Qs-N)>no*SQoJ zA!AfgV|O%G#_)dRBu`yVa{8e$QU;1OTs!c-iGLH;XRLj8k@q~hyyvvc7=F&XQ*R7% zlBX^wxzq9<_4w&bsbas7_dK|~=T2)C{jletdayCbNuIczyk;a%0o=s#e*A*0FkS{!O`NGk! zCQ9pQB3N_1m@)%AkLz*@Imv^|NggwX&lEF*dSiUi=Ii8uCoAKYFU*9m%q_KCS5v(1 zeCn4{r*#J-WQ_S}*&VOz{m2)dxP0MG>n^bRqt^H=KI9}1E+@J3xE=g9(QRP*EI#B5 zPh7rmr>jg;@~G9>3)Y-yU6GSKxSZtCwFE4W1u9@;kS{!S`NGjM-I?H9i7sNFE6YtC zco~?SGNHT`k)fzK~vLV_Sk}o{EeBn9W8nnH~dfs!6PhC#(oURym<#e?Y zy$AWiqsteLhxWi?%I(h2uVZqOr!FTsy2kK(#ru&jJaPHL(Utx}$_>uX??G~s2bYr^ z4`cW}=>3lST}k@Vx1HjLPka>C-`oBy(X*FEK&u8u%otPfusajMRYb_s`;jj^b@{^4 zDkh>JwElZq30P~XKlOU**NdFw(d8sZ*BArA_F3CsIL9DgcyRf`(JDZ<5D&j6N~t$L za*`)5CplWi=n6K!Qm`>fuQ@%v_L75>@3(yp-O*^-lh;6SzE#4xR>6oF!?)M>V>BM7 zbi`Q)C6|40UOeKgCmQ)%O|&ZLM;)aST*+0V?dJwY$QXVOD$zB@2(T$F?LIktxhjV7 z51)#rA0NNBE4ViYuNvkCMv^fq(6T#Y!Sc`y`ZI=n;i=0Pj#fEQMU0D|1Mi_+#sGP3WllLHBcf-zVOuL3r9-~ttGGP_c=MqqsvL2^H{g~z*R68Q~dtlsi7*J)$yk#NA>=9 zSpVmbqjG*4H@}Wy#EjwFuAfD$Ph-I4%&Iq$N&ABO& zFFd+@;dmIsdwB2vo&DRd*TFeI?HlC@-+11;e%2YRC}HnMzVPVsh37mJrE^f@wyuR>YvJn6 z4n&J5+|3T8R&$^zG5i@Nyw0z|BJHYk$?o_w%AYY8KVH+Y@5_Z@#Hq^C>~;OQE393s za~H7uyCbaqp>x;!ZE<>0`uK<&L+Ncg(e&J|UAmCs*YQ;9Rj1tyBW8@LupGNHsXf>mwyWUUFaLcKyqreCmL(gz_$Lroey(zsj{e*bh zf=!Yg&tDVPf2Q(&^YK`NzFG_;%^361vO8W^KPRjtx`sk7@91)Qqh)vQ1g|AxIz`Wk z>y%ibxz8g16;mvw2a}`m9VWha(O2%mp2{@Xs`B1HMm}wp6g81d^dy1 z7vAXdg`;JRkvUu|2G|(n@=jbXZ?sIe2E49?Xioz#;#xu>UwCl&!qGB@*PTYaF}5CZ zPqImyDY5dSfsN~Hd5_EScm{oap1=qh<8rj@P7OHUD(68i?=qLm+i6`5UW#TV<;on^ zIYYki;PQo|WgauYO+-|S0Tx?~JFB)$fBovH*sjlD-ub=e_x21tj141Z46p0=yx*JT z@=jbXZ?tR_y{`8oUwCl&!fCHDG#?9Ssrupl$mN~5T;6ynN^|SilAm+(g$I`}Jg2*k z;XM4D_ZqQPa^npb#>>uJ%xO6?=do_}fo=8q&s`WHW9SUew+d_wzmCb}9b7JNbZ65C zo#EABTb2IvihSXT%NLHWF)G0F6I(wUti$ZI{ch`+F6(_v+^t8~uzvP$zjnPJjF>U1 z(Xy8Kwd>Cka(O2%mv@@2MDH<+dOLl{7am-`a6F!&?b-@uGznf3sj{ltS`_Au27%^j9?L7Rx^Pej} zJ{HQ_zH(LB;fGyBwDw>fUi$sc5O4vwFsz?h7e-^oI0BX!1w;)28$)~0l-BONS^8kb zy=5OBzHTL;y8=cPC#Uxl)+j8*8%E3+(y}{VcLDXLM4s)`<=J)~b>Mas`8~q5fHEH+ za(PFW%R6RT{&#PWP;U(KY^N^IHac2ATG?K=mii8q7V6369bGPOw2Yz8aN#yY)`E>e zp6%fBY;!js1-2+XqpX45Od*$d;&OSTWeh!+)Dp3pB5VxuYzLQT8?Ag4*rLXOwF^|p z<(;@(-e?(P99ZY!GRhdRG03wWU7l^UYzB`2+luiXHh&LFShOg literal 0 HcmV?d00001 diff --git a/examples/pool/main.go b/examples/pool/main.go index 98af7a6d2..56dfcae9d 100644 --- a/examples/pool/main.go +++ b/examples/pool/main.go @@ -9,6 +9,7 @@ Pool Model package main import ( + "github.com/deadsy/sdfx/obj" "github.com/deadsy/sdfx/render" "github.com/deadsy/sdfx/render/dc" "github.com/deadsy/sdfx/sdf" @@ -58,6 +59,24 @@ func main() { } render.ToSTL(pool, 300, "pool1.stl", &render.MarchingCubesOctree{}) render.ToSTL(pool, 15, "pool2.stl", dc.NewDualContouringDefault()) + + // Test ImportSTL + triChan := make(chan *render.Triangle3) + go func() { + dc.NewDualContouringDefault().Render(pool, 15, triChan) + close(triChan) + }() + poolRebuilt := obj.ImportTriMesh(triChan, 0.2, 1) + render.ToSTL(poolRebuilt, 15, "pool3.stl", dc.NewDualContouringDefault()) + + //// Test ImportSTL 2 + //triChan2 := make(chan *render.Triangle3) + //go func() { + // (&render.MarchingCubesOctree{}).Render(pool, 64, triChan2) + // close(triChan2) + //}() + //poolRebuilt2 := obj.ImportTriMesh(triChan2, 0.2, 1) + //render.ToSTL(poolRebuilt2, 64, "pool4.stl", &render.MarchingCubesOctree{}) } //----------------------------------------------------------------------------- diff --git a/go.mod b/go.mod index 522aa671c..f63ceaaa4 100644 --- a/go.mod +++ b/go.mod @@ -5,6 +5,8 @@ go 1.13 require ( github.com/ajstarks/svgo v0.0.0-20200725142600-7a3c8b57fecb github.com/golang/freetype v0.0.0-20170609003504-e2365dfdc4a0 + github.com/hschendel/stl v1.0.4 + github.com/kyroy/kdtree v0.0.0-20200419114247-70830f883f1d github.com/llgcode/draw2d v0.0.0-20200930101115-bfaf5d914d1e github.com/stretchr/testify v1.7.0 github.com/yofu/dxf v0.0.0-20190710012328-5a6d1e83f16c diff --git a/go.sum b/go.sum index bbf40fc22..09a33423c 100644 --- a/go.sum +++ b/go.sum @@ -19,8 +19,16 @@ github.com/go-gl/glfw v0.0.0-20190409004039-e6da0acd62b1/go.mod h1:vR7hzQXu2zJy9 github.com/go-latex/latex v0.0.0-20210118124228-b3d85cf34e07/go.mod h1:CO1AlKB2CSIqUrmQPqA0gdRIlnLEY0gK5JGjh37zN5U= github.com/golang/freetype v0.0.0-20170609003504-e2365dfdc4a0 h1:DACJavvAHhabrF08vX0COfcOBJRhZ8lUbR+ZWIs0Y5g= github.com/golang/freetype v0.0.0-20170609003504-e2365dfdc4a0/go.mod h1:E/TSTwGwJL78qG/PmXZO1EjYhfJinVAhrmmHX6Z8B9k= +github.com/hschendel/stl v1.0.4 h1:DXT5rkiXMUkbKw4Ndi1OYZ/a5SLR35TzxGj46p5Qyf8= +github.com/hschendel/stl v1.0.4/go.mod h1:XQFFLKrq9YTaBpmouDui4JSaxMyAYkpD7elGSSj/y3M= github.com/jung-kurt/gofpdf v1.0.0/go.mod h1:7Id9E/uU8ce6rXgefFLlgrJj/GYY22cpxn+r32jIOes= github.com/jung-kurt/gofpdf v1.0.3-0.20190309125859-24315acbbda5/go.mod h1:7Id9E/uU8ce6rXgefFLlgrJj/GYY22cpxn+r32jIOes= +github.com/jupp0r/go-priority-queue v0.0.0-20160601094913-ab1073853bde h1:+5PMaaQtDUwOcJIUlmX89P0J3iwTvErTmyn5WghzXAQ= +github.com/jupp0r/go-priority-queue v0.0.0-20160601094913-ab1073853bde/go.mod h1:RDgD/dfPmIwFH0qdUOjw71HjtWg56CtyLIoHL+R1wJw= +github.com/kyroy/kdtree v0.0.0-20200419114247-70830f883f1d h1:1n5M/49q9H6QtNJiiVL/W5mqgT1UdlGQ7oLP+DkJ1vs= +github.com/kyroy/kdtree v0.0.0-20200419114247-70830f883f1d/go.mod h1:6oJGQK7VSg3RxSQ7QspgqpCmKjIbAslgT2wBXbFJUZw= +github.com/kyroy/priority-queue v0.0.0-20180327160706-6e21825e7e0c h1:1c7+XOOGQ19cXjZ1Ss/irljQxgPvb+8z+jNEprCXl20= +github.com/kyroy/priority-queue v0.0.0-20180327160706-6e21825e7e0c/go.mod h1:R477L6j2/dUcE0q0aftk0kR5Xt93W7g1066AodcJhEo= github.com/llgcode/draw2d v0.0.0-20200930101115-bfaf5d914d1e h1:YRRazju3DMGuZTSWEj0nE2SCRcK3DW/qdHQ4UQx7sgs= github.com/llgcode/draw2d v0.0.0-20200930101115-bfaf5d914d1e/go.mod h1:mVa0dA29Db2S4LVqDYLlsePDzRJLDfdhVZiI15uY0FA= github.com/llgcode/ps v0.0.0-20150911083025-f1443b32eedb h1:61ndUreYSlWFeCY44JxDDkngVoI7/1MVhEl98Nm0KOk= @@ -34,6 +42,7 @@ github.com/pmezard/go-difflib v1.0.0/go.mod h1:iKH77koFhYxTK1pcRnkKkqfTogsbg7gZN github.com/ruudk/golang-pdf417 v0.0.0-20181029194003-1af4ab5afa58/go.mod h1:6lfFZQK844Gfx8o5WFuvpxWRwnSoipWe/p622j1v06w= github.com/stretchr/objx v0.1.0/go.mod h1:HFkY916IF+rwdDfMAkV7OtwuqBVzrE8GR6GFx+wExME= github.com/stretchr/testify v1.2.2/go.mod h1:a8OnRcib4nhh0OaRAV+Yts87kKdq0PP7pXfy6kDkUVs= +github.com/stretchr/testify v1.4.0/go.mod h1:j7eGeouHqKxXV5pUuKE4zz7dFj8WfuZ+81PSLYec5m4= github.com/stretchr/testify v1.7.0 h1:nwc3DEeHmmLAfoZucVR881uASk0Mfjw8xYJ99tb5CcY= github.com/stretchr/testify v1.7.0/go.mod h1:6Fq8oRcR53rry900zMqJjRRixrwX3KX962/h/Wwjteg= github.com/yofu/dxf v0.0.0-20190710012328-5a6d1e83f16c h1:qgsxLgTXCVH8Dxar36HI5af2ZfinVz5vF8erPpyzM+A= @@ -82,6 +91,7 @@ gonum.org/v1/plot v0.0.0-20190515093506-e2840ee46a6b/go.mod h1:Wt8AAjI+ypCyYX3nZ gonum.org/v1/plot v0.9.0/go.mod h1:3Pcqqmp6RHvJI72kgb8fThyUnav364FOsdDo2aGW5lY= gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405 h1:yhCVgyC4o1eVCa2tZl7eS0r+SDo693bJlVdllGtEeKM= gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0= +gopkg.in/yaml.v2 v2.2.2/go.mod h1:hI93XBmqTisBFMUTm0b8Fm+jr3Dg1NNxqwp+5A1VGuI= gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c h1:dUUwHk2QECo/6vqA44rthZ8ie2QXMNeKRTHCNY2nXvo= gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c/go.mod h1:K4uyk7z7BCEPqu6E+C64Yfv1cQ7kz7rIZviUmN+EgEM= rsc.io/pdf v0.1.1/go.mod h1:n8OzWcQ6Sp37PL01nO98y4iUCRdTGarVfzxY20ICaU4= diff --git a/obj/stl.go b/obj/stl.go new file mode 100644 index 000000000..1cfae6ed4 --- /dev/null +++ b/obj/stl.go @@ -0,0 +1,167 @@ +package obj + +import ( + "github.com/deadsy/sdfx/render" + "github.com/deadsy/sdfx/sdf" + "github.com/hschendel/stl" + "github.com/kyroy/kdtree" + "github.com/kyroy/kdtree/points" + "io" + "math" +) + +// triModelSdf is the SDF that represents an imported triangle-only mesh. +type triModelSdf struct { + // Samples that are on the surface to the owning triangle. + samplesToTriangles *kdtree.KDTree + // bb is the bounding box. + bb sdf.Box3 + // Number of samples to consider + numSamples int +} + +func (m *triModelSdf) Evaluate(p sdf.V3) float64 { + return m.eval(p) +} + +func (m *triModelSdf) eval(p sdf.V3) float64 { + // Find the closest triangle efficiently + closestPoints := m.samplesToTriangles.KNN(points.NewPoint([]float64{p.X, p.Y, p.Z}, nil), m.numSamples) + if len(closestPoints) == 0 { + return 1 // No mesh found: return air + } + // Compute maximum triangle area and + maxTriArea := 1e-12 + for _, closestIthTriangle := range closestPoints { + triElem := closestIthTriangle.(*points.Point) + tri := triElem.Data.(*render.Triangle3) + triArea := tri.V[1].Sub(tri.V[0]).Length() * tri.V[2].Sub(tri.V[0]).Length() / 2 + maxTriArea = math.Max(maxTriArea, triArea) + //if i == 0 { // remove triangles that have >90º when compared to the closest one. + // for j := i + 1; j < len(closestPoints); j++ { + // triElem := closestPoints[j].(*points.Point) + // tri2 := triElem.Data.(*render.Triangle3) + // if tri.Normal().Dot(tri2.Normal()) < 0 { // >90º between planes + // closestPoints = append(closestPoints[:j-off], closestPoints[j+1-off:]...) + // off++ + // } + // } + //} + } + retAccum := -math.MaxFloat64 // 0. + retWeight := 0. + for _, closestIthTriangle := range closestPoints { + triElem := closestIthTriangle.(*points.Point) + tri := triElem.Data.(*render.Triangle3) + triArea := tri.V[1].Sub(tri.V[0]).Length() * tri.V[2].Sub(tri.V[0]).Length() / 2 + // Compute the distance to the triangle (negative for inside the surface) + centroidToTestPoint := p.Sub(sdf.V3{X: triElem.Coordinates[0], Y: triElem.Coordinates[1], Z: triElem.Coordinates[2]}) + signedDistanceToTriPlane := tri.Normal().Dot(centroidToTestPoint) + weight := 1 / (1 + centroidToTestPoint.Length()) // (0, 1) + weight = weight * triArea / maxTriArea // weaker push of smaller faces (0, 1) + weight = math.Pow(weight, 2) // weaker push of further away faces (0, 1) + //log.Println("Weight:", weight, "<--", centroidToTestPoint) + //proj := p.Sub(tri.Normal().MulScalar(signedDistanceToTriPlane)) + //if stlPointInTriangle(proj, tri.V[0], tri.V[1], tri.V[2], 1) { + //retAccum += signedDistanceToTriPlane * weight + //} else { + // retAccum += 1 // Force AIR + //} + //retWeight += weight + retAccum = math.Max(retAccum, signedDistanceToTriPlane) + retWeight = 1 + } + //log.Println("---") + ret := sigmoidScaled(retAccum / retWeight) + //log.Println(ret) + return ret +} + +func sigmoidScaled(x float64) float64 { + return 2/(1+math.Exp(-x)) - 1 // [-1, 1] +} + +func (m *triModelSdf) BoundingBox() sdf.Box3 { + return m.bb +} + +// ImportTriMesh converts a triangle-based mesh into a SDF3 surface. +func ImportTriMesh(tris chan *render.Triangle3, resolution float64, numSamples int) sdf.SDF3 { + m := &triModelSdf{ + samplesToTriangles: kdtree.New(nil), + bb: sdf.Box3{ + Min: sdf.V3{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64}, + Max: sdf.V3{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64}, + }, + numSamples: numSamples, + } + for triangle := range tris { + // Register the triangle's centroid associated to this triangle for later evaluation + triCentroid := triangle.V[0].Add(triangle.V[1]).Add(triangle.V[2]).DivScalar(3) + m.samplesToTriangles.Insert(points.NewPoint([]float64{triCentroid.X, triCentroid.Y, triCentroid.Z}, triangle)) + // Update the bounding box + for _, vertex := range triangle.V { + //vertex2 := vertex.MulScalar(0.99).Add(triCentroid.MulScalar(0.01)) + //m.samplesToTriangles.Insert(points.NewPoint([]float64{vertex2.X, vertex2.Y, vertex2.Z}, triangle)) + m.bb = m.bb.Include(vertex) + } + //// Sample inside the triangle (detail parameter) to avoid problems, + //// e.g. large triangles that are closer than small triangles that have closer vertices + //edge1 := triangle.V[1].Sub(triangle.V[0]) + //edge2 := triangle.V[2].Sub(triangle.V[0]) + //edgeLen := sdf.V2{X: edge1.Length(), Y: edge2.Length()} + //eps := 1e-3 + //uv := sdf.V2{} + //for uv.X = eps; uv.X < edgeLen.X-eps; uv.X += resolution { + // for uv.Y = eps; uv.Y < edgeLen.Y-eps; uv.Y += resolution { + // // Create the "uniform" sample point + // samplePoint := triangle.V[0].Add(edge1.MulScalar(uv.X)).Add(edge2.MulScalar(uv.Y)) + // // Ignore this sample if it lies outside the triangle + // if !stlPointInTriangle(samplePoint, triangle.V[0], triangle.V[1], triangle.V[2], 0) { // TODO: Efficiency? + // continue + // } + // // Add the "uniform" triangle sample to the tree + // m.samplesToTriangles.Insert(points.NewPoint([]float64{samplePoint.X, samplePoint.Y, samplePoint.Z}, triangle)) + // } + //} + } + if !m.bb.Contains(m.bb.Min) { // Return a valid bounding box if no vertices are found in the mesh + m.bb = sdf.Box3{} // Empty box centered at {0,0,0} + } + m.bb = m.bb.ScaleAboutCenter(1 + 1e-12) // Avoids missing faces due to inaccurate math operations. + return m +} + +func stlSameSide(p1, p2, a, b sdf.V3, tol float64) bool { + cp1 := b.Sub(a).Cross(p1.Sub(a)) + cp2 := b.Sub(a).Cross(p2.Sub(a)) + return cp1.Dot(cp2) >= -tol +} + +func stlPointInTriangle(p, a, b, c sdf.V3, tol float64) bool { + return stlSameSide(p, a, b, c, tol) && stlSameSide(p, b, a, c, tol) && stlSameSide(p, c, a, b, tol) +} + +// ImportSTL converts an STL model into a SDF3 surface. +func ImportSTL(reader io.ReadSeeker, resolution float64, numSamples int) (sdf.SDF3, error) { + mesh, err := stl.ReadAll(reader) + if err != nil { + return nil, err + } + triChan := make(chan *render.Triangle3, 64) // Buffer some triangles and send in batches if scheduler prefers it + go func() { + for _, triangle := range mesh.Triangles { + tri := &render.Triangle3{} + for i, vertex := range triangle.Vertices { + tri.V[i] = stlToV3(vertex) + } + triChan <- tri + } + close(triChan) + }() + return ImportTriMesh(triChan, resolution, numSamples), nil +} + +func stlToV3(v stl.Vec3) sdf.V3 { + return sdf.V3{X: float64(v[0]), Y: float64(v[1]), Z: float64(v[2])} +} diff --git a/sdf/box.go b/sdf/box.go index 87503bb13..a59310408 100644 --- a/sdf/box.go +++ b/sdf/box.go @@ -123,6 +123,48 @@ func (a Box2) Enlarge(v V2) Box2 { //----------------------------------------------------------------------------- +// Include makes sure that the box includes the given point by extending it if necessary, returning the new Box3 +func (a Box3) Include(vertex V3) Box3 { + if vertex.X < a.Min.X { + a.Min.X = vertex.X + } + if vertex.Y < a.Min.Y { + a.Min.Y = vertex.Y + } + if vertex.Z < a.Min.Z { + a.Min.Z = vertex.Z + } + if vertex.X > a.Max.X { + a.Max.X = vertex.X + } + if vertex.Y > a.Max.Y { + a.Max.Y = vertex.Y + } + if vertex.Z > a.Max.Z { + a.Max.Z = vertex.Z + } + return a // It is a copy (not passed by reference) +} + +// Include makes sure that the box includes the given point by extending it if necessary, returning the new Box2 +func (a Box2) Include(vertex V2) Box2 { + if vertex.X < a.Min.X { + a.Min.X = vertex.X + } + if vertex.Y < a.Min.Y { + a.Min.Y = vertex.Y + } + if vertex.X > a.Max.X { + a.Max.X = vertex.X + } + if vertex.Y > a.Max.Y { + a.Max.Y = vertex.Y + } + return a // It is a copy (not passed by reference) +} + +//----------------------------------------------------------------------------- + // Contains checks if the 3d box contains the given vector (considering bounds as inside). func (a Box3) Contains(v V3) bool { return a.Min.X <= v.X && a.Min.Y <= v.Y && a.Min.Z <= v.Z && From 140bfa082d8a15390581b8c19b8eeedd4a72d236 Mon Sep 17 00:00:00 2001 From: Yeicor Date: Sat, 1 Jan 2022 10:42:15 +0100 Subject: [PATCH 2/4] Import triangle meshes and voxel-based cache/smoothing --- examples/monkey_hat/main.go | 31 +++- examples/monkey_hat/monkey.stl | Bin 48434 -> 18384 bytes examples/pool/main.go | 19 --- go.mod | 1 - go.sum | 8 -- obj/stl.go | 253 +++++++++++++++++---------------- sdf/voxel.go | 89 ++++++++++++ 7 files changed, 245 insertions(+), 156 deletions(-) create mode 100644 sdf/voxel.go diff --git a/examples/monkey_hat/main.go b/examples/monkey_hat/main.go index 7cbeb1c96..590b9f6d1 100644 --- a/examples/monkey_hat/main.go +++ b/examples/monkey_hat/main.go @@ -12,11 +12,12 @@ import ( "github.com/deadsy/sdfx/obj" "github.com/deadsy/sdfx/render" "github.com/deadsy/sdfx/sdf" + "log" "os" + "time" ) -func main() { - +func monkeyWithHat() sdf.SDF3 { // MONKEY // - Open the STL file file, err := os.OpenFile("monkey.stl", os.O_RDONLY, 0400) @@ -26,11 +27,12 @@ func main() { panic(err) } } - // - Create the SDF from the mesh - monkeyImported, err := obj.ImportSTL(file, 0.2, 5) + // - Create the SDF from the mesh (a modified Suzanne from Blender with 366 faces) + monkeyImported, err := obj.ImportSTL(file) if err != nil { panic(err) } + //monkeyImported = sdf.NewVoxelSDF(monkeyImported, 64) // HAT hatHeight := 0.5 @@ -49,6 +51,23 @@ func main() { fullHat = sdf.Transform3D(fullHat, sdf.Translate3d(sdf.V3{Y: 0.15, Z: 1})) monkeyHat := sdf.Union3D(monkeyImported, fullHat) - render.ToSTL(monkeyHat, 256, "monkey-out.stl", &render.MarchingCubesOctree{}) - //render.ToSTL(monkeyHat, 16, "monkey-out.stl", dc.NewDualContouringDefault()) + // - Cache the mesh full SDF3 hierarchy for faster evaluation (at the cost of initialization time and memory). + // It also smooths the mesh a little using trilinear interpolation. + // It is actually slower for this mesh (unless meshCells <<< renderer's meshCells), but should be faster for + // more complex meshes (with more triangles) or SDF3 hierarchies that take longer to evaluate. + monkeyHat = sdf.NewVoxelSDF(monkeyHat, 64, nil) // Use 32 for harder smoothing demo + + return monkeyHat +} + +func main() { + startTime := time.Now() + monkeyHat := monkeyWithHat() + + render.ToSTL(monkeyHat, 128, "monkey-out.stl", &render.MarchingCubesUniform{}) + + // Dual Contouring is very sensitive to noise (produced when close to shared triangle vertices) + //render.ToSTL(monkeyHat, 64, "monkey-out.stl", dc.NewDualContouringDefault()) + + log.Println("Monkey + hat rendered in", time.Since(startTime)) } diff --git a/examples/monkey_hat/monkey.stl b/examples/monkey_hat/monkey.stl index b18dad2062041560578dec6b9b6ded96bad9bd41..eca6fecae15cd6addef1a20c85a936229bb2ab6b 100644 GIT binary patch literal 18384 zcmbVUd3a6N*WPNL#}Y%$L(F4?$lY5_(FjE;F{H*2Gl?mQp=PC&8bYa52T^Kh&?0v) zEsCO1YAkIGF;lb@9enHDn{&>){rvUa=ebW#a@M=n+WYMB-S0{338O}ii5)ngeAJkc z!^^iAHgLp%fn#bntW&Q}J@NnfjL0Iy&cVmLO`DX^v1j9z--je2&PyTEJwlY>>)U)S z#D>Y`=$*g}z7DJyNHzE?zVZ6*2w@|RF`~NM9~yHXl1kMQ-a# z4Sy|8Y1__l4A;tx9Uk9{2DDtD`=2>08(nHgr%RmFxjU6GGK@Udi@MfG5AJ<-F?(8; z7Dcy9pVD)uUGiYQXl2I8Hhd&)?eKT%srlPD*AUM2eSYrEy+%AE#^RK56s-<<8@HP* z#f1m0m;( zaW!dyyg4P1n&)e%RwvD%yvx(%+1Gz<_V{KZHSe1)??wISL9`HwBMx{*9Gpx?dff^( z^Zk^cjWt5_@gQ1=>~#z1<;COaS-E`rzi(Gk-#&-ssg+OVlyRG=(yN!{(fwZfxNEMcKly%aD}0jH=9NtHlVe z&n#8H*L7xYkauqNq#d#6^vofZWYz_5G1|3yTJWSU@?Ig8kQ~L&qpdeL?YR=F!WXVJ ztT2X&7@3wpGdo}Kyb^lT^J;aS*a^LrEYtfAvF5%zp%S=k6TLv9!EeXnmE ztfLZVm1TrAM>XP~>70%AvHe=t)~s2Gzugwpq>UAVRo|Ziq?xaa#%IIryVErt`mnp2 zw|)Wpcdab1LUduD?|aNA=*_=e=I;`#lB)%C3%JeMU{-#FA8lvF<&bVUtb2Cu&!cTb zdhAPT+|G__@2>fp^Hp|x>st!q-$IF))Kbii}9i!~a;HQ--M@aivklkE0isQz=pPzrh{uR1aN zijK`1;5GdsV)LPrYD(+{Z^(f#1A!Tu=(pPKl6f+|(&2e})8!%&^6CBT6g#x3E;%<% zeiT-VTCXeTT~qrz&+V=4sAly?+`gY|$maD#^`-i?ALzZ2ZRGB8EvQqrd=e3NPy8UO zj1Q&~b7yeBj}GghgYL!9%TIerL|D(}_qomJtapaTjB()m|Mbskz3Ec1s*JZjlAlW< zZCbf#AqM5WCwl}IrCgh9at&Iaa7HxGY#9#P>y~5x_y0*d2K$9O7TS-nuIvnj7=Nmw z?EB?->fY~l8T08g&m+;17JWVVf3o7X%dtCh?3QPX(Tvwe>hGq8QlAs~^|6P|sMs6h zd|Ac)`5JvVexqhjUfYup-%g96IidHxpFLY7F)Kuyh%JXo()7${&d*ufEs^y9g*3txHUL;5P_>&A(p=pqX&PSRXu*(*=i-)ZY6OSi_+!Ec_`@PtzHv7bbkl3 z`v=D0{-NzYh}jD9Y*31h zI!w}o?(Nr2v%mD5r>F~qsJZ-J7Wk;1L@WQ;RD`{xYcs2>i^J*s%97ORs{qfl2J>{4 zhgGS1pMGpxF>qY)Cl5|OSj~_^8BDA4P6pv#RE^-av$^L^){vgJ?c@{ z1M55~w;NKu+}1kOMVOv)tOtMe6yDy9mQ@R3o5ovv%d0y^(XffLG?pAoYa-U^Ua8riPN9Mur|jas zEpcS9Gc*ycI;F{iccLk^Z+0&t&slcO&nomWsGokIjp+1M>`?796CRjV;7f$bNAY{h-mYfbFWn` zxz`rvUW?JWPe{9;;|wQ6M)(kQap^RZ?bVXXxA;ciyZN9QZh zgF@e^ri8e1?X1q9C4!dLSgX-W2$vl+5f@jKP*2i>Xh6D_%rt77X*8`|H1C!tCg^w5 z!sv2_#eunZx_;}?yH?3 zht)b{^;%^!sThy>OJ(z=LhQ~yP<{C5j4V;Hu*IZ0cbE%RHrI=ugkWA&ZDwBdXXZr_ z!K+E`J^pnnY*rOL3Gw!ta_Y(Ua-=->E$-E^V(wMh+$(c{%H{xt_%K_b z3P`I-`zCKO5NO4mjB+`^L!)-espWY++HSZ;zs$obo2^AOkGp#($@a|F?q#-?nMP$Z zw1{B#Q@Lzyi@Aw1{N5Pql6=C8R?LMen@KehmG8`xzvPdlQGK#%L@*DlfQNO^%m{}^ z$gp)|savQ=G9#>PMwq!?Wi!oKJ{~(BKaio!A1KrpGsdIePwP(m2dOKC+wu5>zIw6WeH_Eb zAuE?dPA+m-o@ACYfmu$I-@&XZGumc!7L4pjyUM?$d*&!?@rh1rXFk!1VR~wivRKbL zz0ec!VXR{QTiNV6Gqg@7W@~3=+iNk~lbG2?E3PMw;b*zc^0mm`QVmE=@RlfEMKKGk zY*rr87+Z+p2kWXx=8!iqhm2N;Hsdj~EnQ~&)66>Q^IT11hPNAb@Vt~Tg^mQ2^_V@! zMVs+DuCFHd1V>W;`@c)%m(L!8HTA*)bo_`<$`=E(Fu~Sf$zTN3&Xj6?N17>K} zRYX`$wehrvVwRoIILe@v`2?rjo2b@xh>+vsYFa$C!y#geFSvm9 zbjDFz2xgttI%b_mG3)F{SS>u#VxXO|9b*ge3GeRujUA->d&XxXUc@0=>ks=}m3eCE z^3>KHhiUfpCM%xxSC^AL6kS=Z%{=GyXXTH9^AeVccV1gp%Di!lv`H;|tN|Y4~@m+t6lSSx(}WCEh`V8G}A*qF1iYPOpyW zsu01eNolX6F*c8)l^@A_brel(b4%hKTx*#%QtH%xng7vH`Z*}K2iH{3&L5S3R8$mJ zcay|*5Oc4$%gF|35~%aJnjyFz#1)dK$5uJBbqLj1Rfy{^_mSzc!RB%F&B70il_jEa z1t`Rnz--=_3sb1pAH#Oz-e77L(v=D|r;sh1jKA{Q^|$zbBw<}U#A=79!;wP z*Lpc$XOHt+*<1Qyt;WERo4+I^uLR9iHYZxwCgdw29zyXefuXgTM@NeJ#W* z>qAw>_TO~S$P9@!I(1{ZTyrU$sMSC@H^(ixzH<}m^R>0_tbZ;=maUeBcHDTKa@M~f zn@?;)Rqof8h%jsSZroUPX!jxSs^B>qeIWup@v86dBdI?PC`?N~?5r>=H84qrFKbTk z9r;4WhOd<;6Fbnb>B~7Qd1frltrxEAt_d~nAYygSMY8YQ*JxSiU`Dk4eV{5Zt(^Qd z{flvSPS zO6`S=IQ--%*>7q=DqXt?q2Et-@7JYwzD{v5m7j}ob!b1&y@Rn-zupHL5$G4uLKI&= zPv2`&g_eIFq%f1!Z7=Bw0qtq}ns+6lh4_19oPJG@qQpN&YrKO9GsgNn|9U$mH>dkc zs~9n`^&wh_?`AipC#4?AGkdMokoB=R_VdQLVvZ%5I-D11aV)gG`*rtiKE026k>1DL z3A=l0UslS^%H1QI2sc~K@Rbo7v^bi}z~3CcIm6Y)SO|N@VZDz8y$_e(Mp$zq&Q~_# zVDZ^{;}Rtsd&aQebY73xJ9l?BFy|6RNEdM%vU<$d_AKL;&o@VDn4>_wg9!9gx=JDH z->MPO_zQ-bhBJY9*ZJBm6VFi+<|xqkU?v!k=RpngpriNUnZZ#Dq4vkOOC{K)x_Tc_ z3-KHbiFgLoTE-bwUBR@O0Y|H z5U3VFv=F>Y)v!x-lt*aA`$R`y#A=n%d8oqfF0^!Uc0*Hy2vaA+yHpLkRL2*uULEC; zsovpTs)k*vgFq`(4)F=rjKXx|Xa~Bwrkq9<66S6qcwehwU+XBHutw1nv(CvCsKAKY zw6)x~)=t<#U>ifU5dF(dPfcx=mA=$H4K)?|<^7+=iF9_ z@qBw{343So72u#v#ML?{^#|ojP~PnI4aFAL5AQ!5^%iHfzM6(Tj-%2-?T?AzeUumW zQOhP*dz1M%T&ZArv5%TT~&I3pY)P$$K` zc1IoRo6%G!saId1jk0GlYgY5`UWJHRTa|jp2FUnD@&(!`+sZ=tcdtTBo$#0UVMM4L zQ}20g6#BK+4EA1>_q7`Kwd~8z#Imy235bZ?m!z*(Z=pVGXen*28QE4LVI~$LAK#DS z?~%g~bj&9W)vB_&p3`z60zI(`i9m$}&(^gvG&r9m47Sy`F%;4Cbami4jqKCHiz zpuduKIsR8IWFxrb(k;2=%P!${M%&New&O{NqlX{JF_|UzucMs0au6oZz`8yKx<1bN zJX?xX%*tdxSUV^|JLssX5W&66@#H)0lW^Jxx~@cyHHW`R+A5uU@;7?MEV(7|nP>9d zSTgUp2t2`L4a=a@veoJ8)HtfFp%B88QS+RWwI33+AC5wZ*FKISiBB4ptpnkeq6uyxjfc_5R6E*5vJM* zPx9S(tX7{hluGe4zt7JcPx7%$c%tlUiP~<7+~OSF5?bN<+NYq?prC{K!(IEBa|!GE z2y}g1dVeoMeG#gYgy0^f?H+|ztO|6{{^zkx^4P+EW90yWa==j@VVO)%*;fpqnHTrS zEz|QDDkN4BfR4!O>wbNaPtQQ0XK?DuKO*?8u-;f(JLydVC5NNfLIgf1-Va7Lp<6fV z(Vk&(hSJHj!mL=&;C+#v!Kxoy!(ghA)+X-v9_W%&Zl0gu8?S6D3jvxUtc8y|XiyltC*aaBmB zI*D~l8gxsJHj1@Mj;6@+2ezfpQV1CssP? z7b%@kmt%TjHI;smnhO1zSlV-Q)y6BIZ7}@mx9lV+m%3E#FijvD#3+NNtE#rXq=Tmhwe9OI%qX zS_sx`YtU>zr(3GUb7NZ`-K>v2in5ZNKuHb=cic5cYF3hKP?7_}@@uaRt!q0TuQK#w zUS-VTRR*sO932trq|8zD$EvNURtjq?A3Za4OUU&gT8MvZ-qOV@|XB85Pa;vmopby7mG>ctDH7e~v5R(OWM z3}u;f8k9AT{t6XBsAm#_&x$0R6*<}{u|BHtTGbsq9p*w$`Hcc`b_10M|a|84Jgdb_`Q( z|fbV7uA=Cz^5H*$FC zXcDtKTADj22AV5dVR0;2;ItnUO)}YpiPvA%AQI(lPp2Qq7>KT|rb~INeV%dzs>N@kg zx0&BXodjB8JnrYo6;xqnYqv04i}_+XP?03W$a2#|_Ap!f1+%rN(m}r_;(DEv@_lAf z^DvXjY@MSi!VD438kTEPc3^h!2D5{Rz<8(@5`vjK?J{=(Jhuk&N#afB2nI zf_Fw4;pZ$(zpnPj$Mq?T>-&?rK9l7@D~!SKC=$G*=(G}3Q-WxITUNu{vW`*?l{~0O z^1WrM;4M=whsDWQ>K#|HV`_1jNtNfANj>IcQeDjr*0iZY;x~)E@Mf{2twUW7^VATV z-z=8!X0fALDknkS_oRQkJ!gE#rF9j0*OM#=4K`oA{fMJz4fyx{HpM-ll zjk75Cw41pUD6pLvXoYi$5cBu0l64!5q4uAhvD6}tdKl+g)Kc(V`)s)k?;1;I#owu@ zra;vTdJ>{WMv@#l|;_41{Af{TRQi(tv);EDRX3OJ21X^+V z;GKp>2kv_3G?+{)PHsp=zli2Kbk^GZ4sVXOdQGo58k@MxgLhDuftnERi?=@0Z>||l zi=It4R7jW!Dw2elakP{^`80uKo;qqV<>y8cMMe!Dt4ke?v$ae zJ5z<%-sHZ%t?B2oPE0Opl%p-7DhRFcIq`}tTe1x9*^yqK^w5ZbEfM3f;&=N3-F{jb znzqtI=nG@u^XDD^(E$2nT3fnMxq(J2^lKt6?aZON{n>;@F3FIn-odAdiW}Y)weafQ zTi@^uoj#mUbAvUCHO;4FH=EL`$GPacXXiB{P=AAHZY71{>CP{^Q-i&2BwArBF%ka< DR3WKV literal 48434 zcmbWgdz_YY{>MK;qGTjVsF5&|3^|Nyu4^`_oZ6AL9SI|$O=>nt&Lc@nYzR3I!eG!6 zTFrbuW4{hLEFsn*cF;P8K}&YS@A+!(_kF$V*Z%(a-H*q|^}0XT`}KN#&i6Ibc;AOj zubw<*>Z#+3PMb3M%%VXPPn|UG)G7P)>s!*dq~QP0>-kL!3hq61uk^n!mBu9({x%d8 zv?K)ugJMC!bgl(yK|u+)sjQ%&5ZotM4qby@opqq~Z8{2z#I%#(JW#?NPWjK^$s&jz}+%tOAn z+_Mq6?$N>N)fab&g*VlSfv)bz!;g>6lUC{zg_nkYJGG1LA0D1&>*q7>*nB9;_v!D8 zE=r#Mev{a(&#_{J)_GpX0Bb(PP%(N$QWMwF3g>MgBa-dp=c#~-6HDj zOf3FgNqXpre(~WUDN!vHJDt6*A47MxixC_HUF(zg{&0)mryp+dkEG8xkHGL+t9d^p z*Fx{T!}h~sVM+V6(+#(ZfrriDYOpa@g0&XviRpRb$mFf3HcPpt>(A{7} z#{e5+d=5)#+<>8R(CeF}i@QHg?Z!Nch#F53JI6>fhP78w-p_l`^N{N~8rI%o7&cZd zoLkl>DY>zNVdJByIia;rVta39F3HWH7~Ivkx8rJaE?w{E=A6itp_Q>NDu>1P-nmmW z(}+?JG2pBzhW?nkX36Yt@@VFbr;S{r(3acZXmku%rse0^_oM4+7dKMk;J%;ntmfN0KIfrv^L|8^Wel%7 z9zR>*-s6Y~7dKMkb#E^pxW0$!*5l#*h%U<*Uf2A@@E#9;@>1E(@17SYb=psNC^l9$ zK26klKNyV}V?*8i=@T-g`_H^3+2z3M`03SqZQx;ZzOGf=Gd7Ii7--oYoAV9(anMIE zC*Ag(5Gkwx-;dSkZn!_hh7mG`*IkXDG5lElYfby)%kRGmcdkAuxA%H-okY>AGuI0! zrD%8LI-T-r^FJCdS=Bba+WdRo6{EW`hS7!7E(<2p`ku=#lp z7%^jbj|M!9T?e*SJ-=i~I^^o2xZ7iwr`JcGODr)N1mW9cWPRqyPV9(~33ao#`Q zBEoz0PYVjV5v9I$2Rla07;2eCc|Y&*!B_p#Zd-PZ$Dgwlk1=RWv`KgWp)mgP-R@~Y!CGo`w;Q;GVo~1vDB*r>G97e75P1>7I?-#Z#wsa|KbjMy?>cGJ-K(5d$X7?eE2dHAvF@$)Ikk70ED<+m+j zpka(!aDH^bHf}IN#;Ap5QQpscJoR|T#)~ejh;Izj=TS3SwtjS;Yur|I)!3Eh?$^BY zURljM$46pS6C49<4DWq*_r;B`_uM76EFA{J?+LG~{h|}R39q83ZY!HOrXnsqL+6#> zYb)@umU!<2zC5Yy=FiTL7yM}-wE`{Msb04i^?RT=ACIBcD;tV0E{M+_rB6<^JI?@H zOT71jZto^NyPO_3KlEH?tbdxFI9|7idU=%NvE^U;gw2ke9Gi^uZiC@JGrX?mL7#^;JWF1<$vt7u{a7*Qhwz&vzmF^COhdDPl4pReCEokJy4}L& zyNr*Qy!$*lUQ0VwJuXMrd&5XFhL6&ou$Fl5&F1z^lel;6aLDyU`13=$gYmePtKaAI z-YyG6`yUstZ1`CWdTb|-bd4czJ4^0*adO#F&tDWLuDF9}ud66AG=}Xc+W*2Z8Xdz$ zd5?DZ8DqrTZ!ekL?%a5DpEIddd&|)(2-~>r_rT=FJ6IeyW*H?`hlHo7_)G;_^( z>vK@25x?~a`VGkUyS^ZX)1A$GV|b4RC)6b;-tfoRa?Ex_YfjaYm15^~b@J3uud7x+ zjF2(pA$Q%agk|g+urXeoaZl24;vXYdqJ4Y4uJobyv}y?3vJLPOL1i{sS; zj-lOpX&2Seidx=l#i&mU!K$AJ*6W=~7snn0W^c%7%ovI?_MF~VFPS>^;`s5<%ZTv0 ziqa_brTBiRKIb9E2D-*jlz9j<-stG%$4c)#LsbnXWB9S!cR*#j zvdck@dtcK|QAw-&c*dluto8RWb2b<&7%^is;h|4ZMQNX_?P4#a-bVM@{m)9D*tKW! z-m8B~TF)!p=}J^P-s`hn4DSsiW{h^+wWwt|U9r92ly|QxJLRA6C2MA!mhx$nYr7(r zfz1!jI$TKo6IU*4JhuCDNuS#$h=H#9v5w(9w8m7gHe_=lte09cHpk@qSpiaAr7v4?Q#Z zahIPl{5%crab2=-ufcI}4?W5Fc@|?hb+#vreQ1yC!ox5QZDr)fWZj%QMlJPX*U+ke zuHH3%HF%pibgL~>sb)PdFKhi2)eQQ1f{|v7a54M zyq5YTb`c)UVD&>Y2%|A$XiO~1>o!wwj31Ymr+=HVLp=S7v_qs!V{F+H|fGZ}rg*WAm&exy~eCy(<_-Ldhr!RN-(uBVBCR-OSihCFOU z9-X^+Xn)+Wc>WtZp`n*)`F?nB7)i$P{qSt>J^Rv4V&OJp;^T|1rQLdIZ)-p6{7_4D zLg;1B3@}p1Aks(W*^V)QpSrB9e|{9>p-;8x8~72OW{hI)TGWQ;m3cQl+ddum`N8r0 zzpNnI@AC~k>_ok`<^FJ9td2`M-bl~=KKFj5IqzDoi~n#>(x%sxSl0Ddo|iQT(pt%# z(xkP8it^qo_Lv*q?|63Hss3`JH6M8fk?NBel{qY>%SH`M`(D3w{Iz6hQcy6H($S*8 z_5@H1?nHfa;9q76t73g&bC~(j2x@0ilV^a9;rlV6%W>(1F7Je8%kF}q=O{61i1z*P z-j_apYV32@3*nl5j}n8(JOga~@Z)xN<8f*Ki1xAuFUxIK~F;HHN%}M^lzg+dY1L->C5I>jTpDJv3JO`J;lO zJ~>93F`CfIM}gH6ji~on|G7lEHoC%k($_BbelSuS*VfOdd=wseHrTxRTcU)&5Ne$E z*∨%U0Zm9=3Y)*GK<@fv)*UsTf6lPL)>uW#G{(&x-pUa7%7~*qw6vVdLZdU^HZm z^~WUB)fk#z(G^=Ioi08z4)5_O9$srh59>LM;23U9yr1`Y;-PQCptnZHwZAzPk9JA+ z6K^?KyR?XF2}QlG;tC_l7&>(Qt*(*ZxG1KCw?LOdaFW~NAqHa|b z#iDdR3eV>XBgq(Es~}`-ox8?3@U4BrJNsS`AM4b)VSNv^Mo~>%bv~-K%en6M{>)_f z%V)>UuV2oLd0l;G7{hy9_v-uM&g}=q%l`Be5juBu25Mf)!8Qu=&>Ak}x^m~e(_w$z zK0Y<{%K(FtXMl|{9xT6Rln!6ITe$hie}z-~olK4QQ2WKMA!F_MhYgqB5lKkw1} z!1Ku_T@Q?-KJS8_bPtd&bF&((y|N4s?*}7fjMaGLqrk@Y-jCn7XFBNF9pimJeL=)p zbo=Dm-Z^dkd2<+b-Va917@djkOHtHX6pZb?``2%owtMq;vCpV48d<04?g2w{(FLs2 zXCB!5!H5|{y|<_$u(7@O-ki{+WLroZR-8F7est4i|$O&7c_d6zy8oy-jh9kTL2*7UliC$NP2V$u?~! z#^=X9giakrD^YVihaz3+dp{T+vvZ<-mBs7oE84 zJ?5fdjB2p_7Q=q@q??l~|9pDnN^LzgHV58A^*T$IwK77+@Vb6{gpKj*vZk=~hF3!U zx%#(M@3B0myOgWuBtPeaFDYxOniBtg{hjD|-Q~QOt{C29zdP&7y6@0EPDnNOUaMW& zx`XBG*YRn87@mB0U}L!Upz&g$YwqoMuRSk@_j_|`W7u}s*!bPrYl!yy&U& z4@QzPysqE#dd@aR*Vq3P8g`l%Z#Z-!XNmXdkki#V=nS?~X~>GtlffsP5|ifNvd_J4 zhn%h$-s9$j{v0m+;OJO#Xj^nvbJtE$@3D&d15nWOFpMN)c-_?@W3K`mLenbY|&!M?XUeYj^3hF8=57jPX<83xv;+O|K|hiT^ziwo}CTnu#0eLC0}-tUI0Uz5|X-8WV|xCM-I zw6sc!xO$HY-fLBsaD@>vhSx0*8M^{(jGI2%D&6$y!7;t%A!;>SYOne;jG|V}=USb^ zi@Cyx8Dm(;qWqq(&eeRVZ`s@;4xost5t+7y46<0q#-Va917=BDjavoxP?{m6+5c)j4Lo(>HQ&Q$v zQ;f$7qBOrdgSBq;(kI?BFpMN)G}C(@Wo*wl{g>n7npI1ZJsyehwKmoBwP3ZqJGghc zz8{PRUcOGH%~-hjih&TZlnzo&mO2>9o*!QcTQce`}N8F{^$0 zp96YR>pe6l%Zbu!2-q=V#?XqkDDUSzw%MmK+5h;($zP5-p6H%vZN;^SRoF}mTJUHB zdp{TA0$3bPkZ$;t-^%AA0tBhLa&A3;oui~-yF8v zacX+<)ko_mFJ5^D5qSpK{K{V1x8bsyW7EY46mFoU+0sZ%$<=E%&jVk-|D|E8VH48F ze|rZtdYOAM?~PFiHosXv?VMaSdvwxmlXKvgCars7`Ubk*4@Srs@-TOePe+(yd%xzb zH%W$9O;h1c9Ug7Ds_kA^^)I6^fqEDrV|cAPB8)9;jAJgZPB(w-mSp<2yQ0%TX@|}@ zir1Y-{WGw%%V5Ne;k7h2#-0Z@#)+4mnSNbzZ*pPzKy(&Sw08EVc-yRCjBThLZ{jmE}==-rMD#L{ zHj~KZlxpxYu8S!B8XChHvwEhdT)ht?=i6J#U1RwEb-3<~wEdr}B<^y-tgK8*?M^J^~kcP$HH=awVW%KN_;gVLKj8!{GBzZLaexgP%E z{A7=Je@I5XJ0fy|CfOb7F5>+&MCoOIA9c7dY;otuNzv8%3DRpR${6!uX-(8}eRRxC z$&*DLlF5thcdr&zNrZONQi@)hi&YdDA!BIVeU!018HxsUG? z&id1Lyz?Gayf=pL--~xoNI!b3FkChx65*de1~OxFSrz4!(cmhs)k#4CjF>U>9BNVi zS!6Wz#%SF7#$@Y}Q`5RV`=hZ6jRVkY=Bk&iC1Lp#QcgYZtuTG3lhgaYJDpZ(hBn6V z>+SBcHp$qtPY#X$I5i&VT`|D_;Ygdu^*N#xN{61eJN{cnflNlTn}Cb*7KuQqUvFU zjL|=4QGLP2?g=)=C%avm%sja?E^Dv9j`mtd@xB0E&2uk`R*7a8M#va?Znr3}+lzW* z-1KDE^x~;Ig!^0c`$VrL5B=^!W2pLJ^r4(9jF>Uxk&gnKt~`t}@|1(pE^`JYM}6G~ zzSok+Vsv$~seUC`V+kW>40+_Ez@{q?V_f^s1JgZYN!oOGdogp47BnLu=g{Mwh#7AueoTybkQHqgyGkcZ@c$~5i^Ew@2ZgPNAJB2 zzX_i3MMd(}V~5bz&+vJMaw6D|oA-keGlm~`zgN8X(YN$Zt4lTsHHRJz#iUod$Le3U@-S6ER25?jA}$uZ;J1Y$Jx#^$g#q348YmW4Ag#o$==1 z*>BC1ozaq4J*6wS1SQ`O7)i#c$HStUISGudPfY9i9nW2#+~5AN`01c+(efS}_R0DI zBV-KU5AWf7zVlhThX?yikN?=djCSkQjkat4y`SpM14gnw!|rV8VM^BxI5FLP)Bhx& z{rw5XN&C5pUKZow&t2_#wOHO0cbFWv`01)J{}cU9(GXgcXMiIbI9E z@~{)t`@x7A!=Jlj(9N^K*22FYeSUgpQ+aaX6KnX1URsOM*I7G)QUR``XbkO4{QiTB z;(p(KIqhZ-N^{|8uKw`B^!&4? zg)`gkhKJYHXxbcj4;V3Hc&$o2jN!e%+T!x$@kyiOYn7i6;rmg6M*+Si)XPsT98ZA} zGKO^RP6gN)zE7jB?2?}L=QKU1ME`GwJ||WrM6!;%gS${RT)VrsZyUSpdw9C=k9A_u ziaY~sj5XZr0K18QB`^AIOHy#2`#$+~JpYp8!%;i8Nq3#CpVc(qc}7k_ zqu@QhnYn-Z<$G5p50&mk#6tA7$GgJxe)Flfop#@xfpM?r{+w+2r(MJ#GS2{;pZ8dB z%`eHm$E0bOhdw02kIx!(ZA85{jF2(Bhu8JvcJ_Z7k`E?~jc=Fe?_S5zhib0NC^NwG zxGtw$3a;h4nsU-E_l74I{5E#%^AqpJ(U&{}Yz&mO>CoFA`0eG}%^dpwlW72__jMnnBtFn?q+=8|(`r$G;Le%ys_{tj64uC?p^?vK}p zui8zGHAmhehSvfc zg)6?6&W2;vfQMe*8%B~b)OL&V*ST%J!07jI})KiDE);LuRDW!+rj5O-5ly_`^RurTQMkk z26!nhd2h0Bc{=8& z^OJqQ-QZadGFU6PfT`co*Ym6O8;l^5VQ)o76yKo<#hepePa6Q zX}>ooC2t+1IUkL#J@5K{jG?}qJ}p2EM$8!f(Y3aZ1{-?}*cd~e+pX;2E34xJZ61d2 zpFby}TSc8l)t~t=8Z(AJ!$+W7L3EuM;A*bRD81KSnSB5Gf>6K3cr^WfuHgL~B5D{> zKW;EW#;8EcqWt*yar^D)eRVsNs5D*N2j#HX8aHG7?f$LP|D4z(a-~*38&smD zQT2YKs28Jza_3E_#iy6-7*4)>S!n&)AkWC@%ER~X>B~E%|LQgX-G{DGt$ldRG@Aq*$Hofbj(&W$~?b7vsR`Y&-O~Z&8!+VUuubgqxF7oSm z#X%F(-W~5RTlr!Wt@8WMkJ}jhjs4%&V`T4kG*u`n(<@3rY z>;vynimv8ZCyubKP#7^|cr9ULc<)zdl%zM;UKj3N_b=6c4m9WLbv0!wcp+DfI*gbx zH0~DV=VL1M#u#2WEM587qQ=vA=}xU*qu#@>G5vq9U-AFHw*KF;e*JjA%AB|Qb=GtJ z(!q~?)5vuLkHKkc3$@xc|7vZH5i`c%oR*@zpZC~u&Cg|READUDa?OS4X#Lpk^d5dc zzRUj;{_dO~%D9SwN1g#zOXR1sOSl$%{H#Bw$1HrSv23$DsZo2y(BBB_e6%R3 z<6R%`61vp>MGUm^45zDeUTi6SwoQ-p!0neee)NVuZM=uh#C1`utI>vA88Kt%+_fn0 z=RNFT>+*C@^Bim~yRhB~Ri2pPj` z-GG9z=Yx%L&d*20Yx~@uJTj#%wYQ_S5}g`~*S(YawbV_gz(_NO*SZ}AW8VojM&Y5q zjsJbIUE2JF_})X$qU&11>XV)PeYt7>Q+okFJc`0KF}GZb<|A#NU&C- zcJ2M&?H8Wi;eyz1pYH4yzYla?8KW9p52FUg9-Cf}oZI`H_`B&3^R5~#yQ7n05m7p$ zyx+}_yb=EI-SKhQaeIm3w5EZLQ3GDrKiC&x1AX0{PHaCH0QQEy{`9zF)(BJ_WJ$k{my89I(fZmVtk?J z_KoY0l|M@~ADVNGRXF;cF#f*NgXSbclSm~_IGn}7aW378kt3WTmMk~{dJvfJ(P;Gq~!+VTGMSfc0 z{_}D2i9dy3h8M-Q@9l)nSjudobe{XYQjzmeefN9kC*AMu8t;8I5hFSVSRVRZ@@u+% z<=$z(S*4M}-C3~gH1K=Hd&5w4j^VXt;bCm={cU}7GNZIAzIE-JeG0=U6)@c`c4~?bm!DoAPPW!GsIv#OVI~Z5f53O*mYOm|h_Uhlw55x93 zAzt#!t74#)XMoqh)K3uJ`=1>O(tgvAiJNbC4R`%#hWGOxFk;5=9{#h$pY2Uo&rd$x zYEb<3iq7=G??J!syf=)HF}$wdcYZ&9G(42`8~JqErQt$)Dvvt)qu*H_k;4Uyi0apB z6sU*Mm@$rkCC0j_92Q$1#yEAemZase?ZWdroJ8${H2Ync*R7*oKRZ-0(lA2C@LCH{ zSb?9S>cGZ$cKD=p#__Yls9k2zUhnY;ngx{ku>9}(X6|)n?ET#};r)M}FNV{71Z)iN z@vp}!)3Gm}RW`YL3-r8JE$=%}TJX^4yjFf2uKU&=8?U_clmY!a6!Omhp0YM(4DUB$ zPEBZd@c87nck5qYz1F3?*Gg=r2-k95O)?@Z3QeWDmKbJz2gF}&a3-up4! z`_p!1zr1=946ikY_u5@5G|W* zulF9jWc~HyJ$#>5ZP_8+?yynmYww=H>E-9#d-(M>zIT`K#laV*e=qLMN#?cuTJmcQ zMv^i7-1@eA?^V5DO$N4^kTx{EgNL7IzplI=jF2(Bho5IZAFn+)Jl;2K){=&8cVKi! zP%0R6?J7T3L#Wr>&!*gXFaLsd$Zr}?d~zQ#qGNz!utQ7 z>Gydr)U>+_m;o3eWB73&g08is7uXnm9+=bko3rmL8~mq>sO^eYv9$2>9_mK{^@Y^K zD9ae$qbmx=mY*@w_xDc!_Q^r3Ok>r|6ffzw+}3N6o8Gb{oa{+okTEV#^^w;-SK{f)SJ>5Yu`y8 zzWco};rk=0^}2eZDODZ1!rj4oriKwR#v)o~cZ7}M{T}aqUvlMk!{euSD!Z5;AL$N5 zw;fmS2P0$*uPbZ}KW?x5G_&y^laGmO59@@U*VQxVQPgP`D=3NjL7!lhWsC}R?amyq zF}&X{edmOW&-*bv_5R45hxOL$7U8i6_cb3dl8jM=uHEr|@-QX7bs@TMT|i5%(wx_# zSDeG%gKu_7-8Va+Weo4(y|=sY>(FcJ331ws3DJIhOyY85-7~U_J{oB5EfxBvz7&bmyYku7L zW{1>$vqQ)jzE8D{*+}!P3&DNs0(zRM&$uq4bO&#iuJ6Y;J0$L#9q=%Q?~~1kTFAF9 zG`eqHphhcvF;{uy*W?O$qE8+5iEnm@?wcJ##%KXs3s-}U;d1uV zR*idp*1*HjZQ|?#^75QqWjhbBJIQ(50;0<$NQ1jFm>PTfR-`l9a=4hQc@%+Rq#47sbjLiz@aL8LP7B&>W2GpYn|3*@p0{w`yRo9i#zc2) zti+fTT3ah6_uU*~P@@r5R5L|uOi?!HYOR}d?yd~24BfFg7sLC#`G=9|n}1lIbm*hM zq4&BiMCMl_*j6Hpm@!tvvO6tcV|YKZ%9pvU@@Vy>b(1K1_2N1uhqY@jppdyfxXktF z8lx24k%;LVFvu#OxUBLyj|oJr8%eJ6E2ZApWUdb`bA8T3j5=^zu4`$oy>x?>p;{9O6jNx@_&^3nU+?2k!c6hpfpDwY>K3B2&YI7c%u}bhk=&uBO zKNvA%c-=}OjZq7hhm?-_X;SjW@!Q88cOObkMbdhnmvg5Xtlk!aUk2+8F6UDkM#va1 zbJy-j&lnY8eKJUitn!J=DxdP#RNN^-t2?*|tmiiGLFW44GS?Gf4DaE+pZ{{-bm*r$ z#6xd-3&)=58e=3_H$Bm${yZt`sr6hxaC{eCo2wqpR^*1EU12$mJXOk-0v&%=JVVLq9v4pPpi@ zg=CdaT~>K%bGKXCI_ElrJAl0(nd_s=Tpu!q)|B~4#aejB)w`sWck{1 z4#J}rJPO`+bXnzDr{(CXKg+@4F&s(_qw);Mjxjwqg^|aR*#bAAwXeD}&zP~*wxxCNu z@r8H3A)v}&j*O$4>^=KJm z1-O(}HKXGF$SU9HvdX8~s#av44K_BJ>&sl`dNgexbOZOMG!gCn$SOauILj)JhwVu2bZ~?Xk%!PE2K&!X<==GnP= zvB_K?T;}?mt{6?=1LVPP^t>NgXaarZjuoFOR zJ4*`i8H~rRT*+J?UFLdvYz%n}Axita0Bnqb!)GNw4LmLWx6@ZMgJR_kTANxo z`*R&j(Mz8)<8rtc*C`Yj12cxA?2c#X3^t`F_Mesv>V0%PX#WoI^;3jCG4+|}{WSkQ zsnfGBjF2%j2X@Enw%}n(`Yt_}Ro;12=X58aI~%OCL~DS|^}%JXcUsM0eb(w^n~k?I zdYv^lxnl0w@wNSD@w}_Cm-|X0YVp%Z3#)%(YwTczjG^}0om#LlDq)$@AFjGDId^CB zk=@dX+UA_rx|SqaOKQN<*FJ<1GKM}U?2d1{?~~?f<5}f($6Dq4MCST|O&iTzzmR@t zpU=bB_sMw>;kCSnutuaWh0OH>T;_T=K7QOhL)c}NcVp$}U_%D04mZ6v922!fJ4r8V z6>Bh;}7H^DnlT-Va917+zO>GR88nDUnsa%w?6&>8cgpL(e|OAai}A%Uqw+ z6+@x#v&t7aKhGvdeByG%=iIfH{D_L7J}sc0JoLfkp{G@9i5Pkb*X6J= z$Pu5o9Pzo9Xr^kx!(r6q7{XE_4}EZX=<$%27NDRQ0X%T%Kc$ zgssn~a{Ma5nps6=jB-5ej@ZWVN)vzGGCAdrPO;~&!}0JV>UEoP9z}S-2pK~wQ69SE z{fekpp9-iaM||pX#5=8Kuv#^qD8CZPLmyopdgsxc(_OcAVcE)mvElV({?@mLx|zqL z=e4{?Eq*pPFha(tMc3}gBi}0LL5}z`mm?l6i}K@BPQ8sfdFUHm9(t$i=cAl@V`z`x z^1)R0`eiu#C}CR}z8^3`#_;{{46pQmgD+_ubLO_OaT|SN>6PDiv|U(hsUT@>mG&Zx zGRNSq+N(Q$r+TG-&lr$?e1G>ix@;~w{sfS&ol(VLeMU`35k|}y#njjxts!H0KXSw; zE=PQ>?b>UlNh`Y+>^;asA6*`L+MnO&&cl0?BR+9C;_*;RRub($sl6Y0=!45ckB2c9 zI=@mB)IYV59Pz2k5uddFhE8|9uJx5i^D#_X4zS zjp>vt#Jij;IpULz=ZM$p(-V!JtF<5XT&me64}I8p9(wI_eIkq11UAN_75&1$ubvjq z?{^H2`qbAbj-oV!^9-ds-cS7yL;Dm) zS;mlu-SN6L)SJ?#t9pdZKbjmjU3EvI(aD~?yfstg4aDP;hF7}k;tugzzN@g4B_kB%Spc!CJ8Ta;$M;&H-Da+DQ?wFa=IZ^(4If;N_>eKWIln@%wQ$wv z4nYl`x(_2{jQWt>sRJ8h1=y56zOQq#dCv*)zVfHx zd)*q|OSgsV9E$p+zQG6?V@}BK)PRlQ{l5Rh4+A^?cY6HO?i1mA-D2KLcRANeitSn$ z4H=^{WOs_e#_)dRZ%}yhlB1rt}Oy8@T-K&O?4Wsh4rpcnL4%N{;x%<%kcN z*4tpUO1)hUHU|0IgUjE}-Fy_-qS~dcyIH>Eh;MW`;?XjOU$5=((*76Y`?nqpn;i6q zc<=?vu;KT!-xJ;gMv^hSmiO>`ko@hzx9!>bjV??FLJ~uE=N2O#_(rX z2kONZ=f;cuAO$?1IECN;kCTS>X0dszrE4rZ+9N6zz5J0 zJ(KD*(i6u*FgfA}x*YLNOXs0o30?~L!`(aSTyj`2S-zI%^yJha}5V1D22i1e<8wz2pTeZPTbNV*!E(j1UxVRuTv z#_)bq&p15YzolDz;*%d4MXy_qp7pi@EG6BC5i>>wx^~BpkN10`v1c6n^IwxWC-hDA z{Mb5!E6~!DhS+Ptn#G;LtH3bQjL}T*?M@BY7+z_Q=Z=fv!KahQK2?TTudCV89F*d% zvFZ-)4Th0s42_lDDFPeA`xX6M6S8xf~dypeOaXI4AHHP=_-sEpj zUH*1-bux5AcMZ5BSZ9j&BS(C6IpXm!#scS854IMPzddpJ+tDjRSNfWtwM2df_I~7u z4=zVM(Z*1;`KH4pquX#I@~jF>UzQEPWJ z$HrI$HlU@ZeJ*WCn(xqj)KWBxYX71f*4)&9O$kQK7&S!N9ckMf)aJ01$lsp2{Ot+tq;@W) z65PUdJf#|})htEuJT=zC3UDD&n@t$!i35WE$*t5kUP0S^O{BR;qs z@zfYYQTmjdTPBn-$lsp0{OwM6K6o}oqpa5=@H1d?#0Qrno|=ud#50~Mnw4B$@>AK| zE@#JeQh6=!p}8qWO{)|}$QYH>*q!;{jkW5#@i!$8O>3VlXt1pEnz2PfMBc7U>l;xCa$})=9aTUcFRd4T`zCUtp6-BP&3K;DyP7R8AH#ZcE>-X458kX$lsp0{OxETfS%TBJ!NAn zh8*$1<%ka%L*p|Vb|FPCTVv#JFLU|Z(N&KX`59QNa8iy@3%&qMj`-km#1m}{>CWT* zX&W$VrvF&>@8gChQ%W|L@!dS+_YCDk*f1J1#%E~S9bse0&sufjGquTQdksqO>DMXe z-hr!H;&oS{Z(|4}WQru1FuA!T1Z)-67E_mSx6S#&O1^6fd6*LU+f$dn9gQxOI`rhXlB?#(dypeOx*YK_ zV`$%LXU(Q~?@rI&m6U$`oAlOAwXXd7`HV;u#DPVX0v_)=ryb%kvQ%fn_6M$8!UushN< zhWESa(314E`*w-Dch`9HHM<}>|>{rF|I7_ zA70t)ytHf07wj?ryO&{jbcdnyy&l{ZTmpuXWQ=+|?9MQ-F`C)QJHyZjEbVf6sGpZ6 z2bWEO?>*LVH{WxxwH8K_G5nbLvGN||R!?1Sb$ps=O);9I!CEoZV1IUzzdgGA?X=Yx z+IR9?cZ!15N8z(yIX9g$x-@jz_7}#=>yE*roVo=>X`dIs(3uD$W(?`t9sQoc7^A_~ z!dYFWrM<42l-#qV3tC=x1bP+J)#W_w`F`7F2{F&(e z$lo4a{&sYY;m<4YM{f1h~3f+ zhozTJ=@xJkYs>coM$8z#z0E`$qZ?|bbkcseC7sSZGS1qrEwx@(<0Qs-N)>no*SQoJ zA!AfgV|O%G#_)dRBu`yVa{8e$QU;1OTs!c-iGLH;XRLj8k@q~hyyvvc7=F&XQ*R7% zlBX^wxzq9<_4w&bsbas7_dK|~=T2)C{jletdayCbNuIczyk;a%0o=s#e*A*0FkS{!O`NGk! zCQ9pQB3N_1m@)%AkLz*@Imv^|NggwX&lEF*dSiUi=Ii8uCoAKYFU*9m%q_KCS5v(1 zeCn4{r*#J-WQ_S}*&VOz{m2)dxP0MG>n^bRqt^H=KI9}1E+@J3xE=g9(QRP*EI#B5 zPh7rmr>jg;@~G9>3)Y-yU6GSKxSZtCwFE4W1u9@;kS{!S`NGjM-I?H9i7sNFE6YtC zco~?SGNHT`k)fzK~vLV_Sk}o{EeBn9W8nnH~dfs!6PhC#(oURym<#e?Y zy$AWiqsteLhxWi?%I(h2uVZqOr!FTsy2kK(#ru&jJaPHL(Utx}$_>uX??G~s2bYr^ z4`cW}=>3lST}k@Vx1HjLPka>C-`oBy(X*FEK&u8u%otPfusajMRYb_s`;jj^b@{^4 zDkh>JwElZq30P~XKlOU**NdFw(d8sZ*BArA_F3CsIL9DgcyRf`(JDZ<5D&j6N~t$L za*`)5CplWi=n6K!Qm`>fuQ@%v_L75>@3(yp-O*^-lh;6SzE#4xR>6oF!?)M>V>BM7 zbi`Q)C6|40UOeKgCmQ)%O|&ZLM;)aST*+0V?dJwY$QXVOD$zB@2(T$F?LIktxhjV7 z51)#rA0NNBE4ViYuNvkCMv^fq(6T#Y!Sc`y`ZI=n;i=0Pj#fEQMU0D|1Mi_+#sGP3WllLHBcf-zVOuL3r9-~ttGGP_c=MqqsvL2^H{g~z*R68Q~dtlsi7*J)$yk#NA>=9 zSpVmbqjG*4H@}Wy#EjwFuAfD$Ph-I4%&Iq$N&ABO& zFFd+@;dmIsdwB2vo&DRd*TFeI?HlC@-+11;e%2YRC}HnMzVPVsh37mJrE^f@wyuR>YvJn6 z4n&J5+|3T8R&$^zG5i@Nyw0z|BJHYk$?o_w%AYY8KVH+Y@5_Z@#Hq^C>~;OQE393s za~H7uyCbaqp>x;!ZE<>0`uK<&L+Ncg(e&J|UAmCs*YQ;9Rj1tyBW8@LupGNHsXf>mwyWUUFaLcKyqreCmL(gz_$Lroey(zsj{e*bh zf=!Yg&tDVPf2Q(&^YK`NzFG_;%^361vO8W^KPRjtx`sk7@91)Qqh)vQ1g|AxIz`Wk z>y%ibxz8g16;mvw2a}`m9VWha(O2%mp2{@Xs`B1HMm}wp6g81d^dy1 z7vAXdg`;JRkvUu|2G|(n@=jbXZ?sIe2E49?Xioz#;#xu>UwCl&!qGB@*PTYaF}5CZ zPqImyDY5dSfsN~Hd5_EScm{oap1=qh<8rj@P7OHUD(68i?=qLm+i6`5UW#TV<;on^ zIYYki;PQo|WgauYO+-|S0Tx?~JFB)$fBovH*sjlD-ub=e_x21tj141Z46p0=yx*JT z@=jbXZ?tR_y{`8oUwCl&!fCHDG#?9Ssrupl$mN~5T;6ynN^|SilAm+(g$I`}Jg2*k z;XM4D_ZqQPa^npb#>>uJ%xO6?=do_}fo=8q&s`WHW9SUew+d_wzmCb}9b7JNbZ65C zo#EABTb2IvihSXT%NLHWF)G0F6I(wUti$ZI{ch`+F6(_v+^t8~uzvP$zjnPJjF>U1 z(Xy8Kwd>Cka(O2%mv@@2MDH<+dOLl{7am-`a6F!&?b-@uGznf3sj{ltS`_Au27%^j9?L7Rx^Pej} zJ{HQ_zH(LB;fGyBwDw>fUi$sc5O4vwFsz?h7e-^oI0BX!1w;)28$)~0l-BONS^8kb zy=5OBzHTL;y8=cPC#Uxl)+j8*8%E3+(y}{VcLDXLM4s)`<=J)~b>Mas`8~q5fHEH+ za(PFW%R6RT{&#PWP;U(KY^N^IHac2ATG?K=mii8q7V6369bGPOw2Yz8aN#yY)`E>e zp6%fBY;!js1-2+XqpX45Od*$d;&OSTWeh!+)Dp3pB5VxuYzLQT8?Ag4*rLXOwF^|p z<(;@(-e?(P99ZY!GRhdRG03wWU7l^UYzB`2+luiXHh&LFShOg diff --git a/examples/pool/main.go b/examples/pool/main.go index 56dfcae9d..98af7a6d2 100644 --- a/examples/pool/main.go +++ b/examples/pool/main.go @@ -9,7 +9,6 @@ Pool Model package main import ( - "github.com/deadsy/sdfx/obj" "github.com/deadsy/sdfx/render" "github.com/deadsy/sdfx/render/dc" "github.com/deadsy/sdfx/sdf" @@ -59,24 +58,6 @@ func main() { } render.ToSTL(pool, 300, "pool1.stl", &render.MarchingCubesOctree{}) render.ToSTL(pool, 15, "pool2.stl", dc.NewDualContouringDefault()) - - // Test ImportSTL - triChan := make(chan *render.Triangle3) - go func() { - dc.NewDualContouringDefault().Render(pool, 15, triChan) - close(triChan) - }() - poolRebuilt := obj.ImportTriMesh(triChan, 0.2, 1) - render.ToSTL(poolRebuilt, 15, "pool3.stl", dc.NewDualContouringDefault()) - - //// Test ImportSTL 2 - //triChan2 := make(chan *render.Triangle3) - //go func() { - // (&render.MarchingCubesOctree{}).Render(pool, 64, triChan2) - // close(triChan2) - //}() - //poolRebuilt2 := obj.ImportTriMesh(triChan2, 0.2, 1) - //render.ToSTL(poolRebuilt2, 64, "pool4.stl", &render.MarchingCubesOctree{}) } //----------------------------------------------------------------------------- diff --git a/go.mod b/go.mod index f63ceaaa4..4520fdd0c 100644 --- a/go.mod +++ b/go.mod @@ -6,7 +6,6 @@ require ( github.com/ajstarks/svgo v0.0.0-20200725142600-7a3c8b57fecb github.com/golang/freetype v0.0.0-20170609003504-e2365dfdc4a0 github.com/hschendel/stl v1.0.4 - github.com/kyroy/kdtree v0.0.0-20200419114247-70830f883f1d github.com/llgcode/draw2d v0.0.0-20200930101115-bfaf5d914d1e github.com/stretchr/testify v1.7.0 github.com/yofu/dxf v0.0.0-20190710012328-5a6d1e83f16c diff --git a/go.sum b/go.sum index 09a33423c..2b9d1924e 100644 --- a/go.sum +++ b/go.sum @@ -23,12 +23,6 @@ github.com/hschendel/stl v1.0.4 h1:DXT5rkiXMUkbKw4Ndi1OYZ/a5SLR35TzxGj46p5Qyf8= github.com/hschendel/stl v1.0.4/go.mod h1:XQFFLKrq9YTaBpmouDui4JSaxMyAYkpD7elGSSj/y3M= github.com/jung-kurt/gofpdf v1.0.0/go.mod h1:7Id9E/uU8ce6rXgefFLlgrJj/GYY22cpxn+r32jIOes= github.com/jung-kurt/gofpdf v1.0.3-0.20190309125859-24315acbbda5/go.mod h1:7Id9E/uU8ce6rXgefFLlgrJj/GYY22cpxn+r32jIOes= -github.com/jupp0r/go-priority-queue v0.0.0-20160601094913-ab1073853bde h1:+5PMaaQtDUwOcJIUlmX89P0J3iwTvErTmyn5WghzXAQ= -github.com/jupp0r/go-priority-queue v0.0.0-20160601094913-ab1073853bde/go.mod h1:RDgD/dfPmIwFH0qdUOjw71HjtWg56CtyLIoHL+R1wJw= -github.com/kyroy/kdtree v0.0.0-20200419114247-70830f883f1d h1:1n5M/49q9H6QtNJiiVL/W5mqgT1UdlGQ7oLP+DkJ1vs= -github.com/kyroy/kdtree v0.0.0-20200419114247-70830f883f1d/go.mod h1:6oJGQK7VSg3RxSQ7QspgqpCmKjIbAslgT2wBXbFJUZw= -github.com/kyroy/priority-queue v0.0.0-20180327160706-6e21825e7e0c h1:1c7+XOOGQ19cXjZ1Ss/irljQxgPvb+8z+jNEprCXl20= -github.com/kyroy/priority-queue v0.0.0-20180327160706-6e21825e7e0c/go.mod h1:R477L6j2/dUcE0q0aftk0kR5Xt93W7g1066AodcJhEo= github.com/llgcode/draw2d v0.0.0-20200930101115-bfaf5d914d1e h1:YRRazju3DMGuZTSWEj0nE2SCRcK3DW/qdHQ4UQx7sgs= github.com/llgcode/draw2d v0.0.0-20200930101115-bfaf5d914d1e/go.mod h1:mVa0dA29Db2S4LVqDYLlsePDzRJLDfdhVZiI15uY0FA= github.com/llgcode/ps v0.0.0-20150911083025-f1443b32eedb h1:61ndUreYSlWFeCY44JxDDkngVoI7/1MVhEl98Nm0KOk= @@ -42,7 +36,6 @@ github.com/pmezard/go-difflib v1.0.0/go.mod h1:iKH77koFhYxTK1pcRnkKkqfTogsbg7gZN github.com/ruudk/golang-pdf417 v0.0.0-20181029194003-1af4ab5afa58/go.mod h1:6lfFZQK844Gfx8o5WFuvpxWRwnSoipWe/p622j1v06w= github.com/stretchr/objx v0.1.0/go.mod h1:HFkY916IF+rwdDfMAkV7OtwuqBVzrE8GR6GFx+wExME= github.com/stretchr/testify v1.2.2/go.mod h1:a8OnRcib4nhh0OaRAV+Yts87kKdq0PP7pXfy6kDkUVs= -github.com/stretchr/testify v1.4.0/go.mod h1:j7eGeouHqKxXV5pUuKE4zz7dFj8WfuZ+81PSLYec5m4= github.com/stretchr/testify v1.7.0 h1:nwc3DEeHmmLAfoZucVR881uASk0Mfjw8xYJ99tb5CcY= github.com/stretchr/testify v1.7.0/go.mod h1:6Fq8oRcR53rry900zMqJjRRixrwX3KX962/h/Wwjteg= github.com/yofu/dxf v0.0.0-20190710012328-5a6d1e83f16c h1:qgsxLgTXCVH8Dxar36HI5af2ZfinVz5vF8erPpyzM+A= @@ -91,7 +84,6 @@ gonum.org/v1/plot v0.0.0-20190515093506-e2840ee46a6b/go.mod h1:Wt8AAjI+ypCyYX3nZ gonum.org/v1/plot v0.9.0/go.mod h1:3Pcqqmp6RHvJI72kgb8fThyUnav364FOsdDo2aGW5lY= gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405 h1:yhCVgyC4o1eVCa2tZl7eS0r+SDo693bJlVdllGtEeKM= gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0= -gopkg.in/yaml.v2 v2.2.2/go.mod h1:hI93XBmqTisBFMUTm0b8Fm+jr3Dg1NNxqwp+5A1VGuI= gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c h1:dUUwHk2QECo/6vqA44rthZ8ie2QXMNeKRTHCNY2nXvo= gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c/go.mod h1:K4uyk7z7BCEPqu6E+C64Yfv1cQ7kz7rIZviUmN+EgEM= rsc.io/pdf v0.1.1/go.mod h1:n8OzWcQ6Sp37PL01nO98y4iUCRdTGarVfzxY20ICaU4= diff --git a/obj/stl.go b/obj/stl.go index 1cfae6ed4..e808bb9db 100644 --- a/obj/stl.go +++ b/obj/stl.go @@ -1,167 +1,176 @@ +//----------------------------------------------------------------------------- +/* + +Closed-surface triangle meshes (and STL files) + +*/ +//----------------------------------------------------------------------------- + package obj import ( "github.com/deadsy/sdfx/render" "github.com/deadsy/sdfx/sdf" "github.com/hschendel/stl" - "github.com/kyroy/kdtree" - "github.com/kyroy/kdtree/points" "io" "math" ) -// triModelSdf is the SDF that represents an imported triangle-only mesh. -type triModelSdf struct { - // Samples that are on the surface to the owning triangle. - samplesToTriangles *kdtree.KDTree - // bb is the bounding box. - bb sdf.Box3 - // Number of samples to consider - numSamples int -} +//----------------------------------------------------------------------------- -func (m *triModelSdf) Evaluate(p sdf.V3) float64 { - return m.eval(p) +type triMeshSdf struct { + tris []*render.Triangle3 + bb sdf.Box3 } -func (m *triModelSdf) eval(p sdf.V3) float64 { - // Find the closest triangle efficiently - closestPoints := m.samplesToTriangles.KNN(points.NewPoint([]float64{p.X, p.Y, p.Z}, nil), m.numSamples) - if len(closestPoints) == 0 { - return 1 // No mesh found: return air - } - // Compute maximum triangle area and - maxTriArea := 1e-12 - for _, closestIthTriangle := range closestPoints { - triElem := closestIthTriangle.(*points.Point) - tri := triElem.Data.(*render.Triangle3) - triArea := tri.V[1].Sub(tri.V[0]).Length() * tri.V[2].Sub(tri.V[0]).Length() / 2 - maxTriArea = math.Max(maxTriArea, triArea) - //if i == 0 { // remove triangles that have >90º when compared to the closest one. - // for j := i + 1; j < len(closestPoints); j++ { - // triElem := closestPoints[j].(*points.Point) - // tri2 := triElem.Data.(*render.Triangle3) - // if tri.Normal().Dot(tri2.Normal()) < 0 { // >90º between planes - // closestPoints = append(closestPoints[:j-off], closestPoints[j+1-off:]...) - // off++ - // } - // } - //} +const stlEpsilon = 1e-12 + +func (t *triMeshSdf) Evaluate(p sdf.V3) float64 { + if !t.bb.Contains(p) { // Fast exit + // Length to surface is at least distance to bounding box + return t.bb.Include(p).Size().Sub(t.bb.Size()).Length() } - retAccum := -math.MaxFloat64 // 0. - retWeight := 0. - for _, closestIthTriangle := range closestPoints { - triElem := closestIthTriangle.(*points.Point) - tri := triElem.Data.(*render.Triangle3) - triArea := tri.V[1].Sub(tri.V[0]).Length() * tri.V[2].Sub(tri.V[0]).Length() / 2 - // Compute the distance to the triangle (negative for inside the surface) - centroidToTestPoint := p.Sub(sdf.V3{X: triElem.Coordinates[0], Y: triElem.Coordinates[1], Z: triElem.Coordinates[2]}) - signedDistanceToTriPlane := tri.Normal().Dot(centroidToTestPoint) - weight := 1 / (1 + centroidToTestPoint.Length()) // (0, 1) - weight = weight * triArea / maxTriArea // weaker push of smaller faces (0, 1) - weight = math.Pow(weight, 2) // weaker push of further away faces (0, 1) - //log.Println("Weight:", weight, "<--", centroidToTestPoint) - //proj := p.Sub(tri.Normal().MulScalar(signedDistanceToTriPlane)) - //if stlPointInTriangle(proj, tri.V[0], tri.V[1], tri.V[2], 1) { - //retAccum += signedDistanceToTriPlane * weight - //} else { - // retAccum += 1 // Force AIR - //} - //retWeight += weight - retAccum = math.Max(retAccum, signedDistanceToTriPlane) - retWeight = 1 + // Check all triangle distances + signedDistanceResult := 1. + closestTriangle := math.MaxFloat64 + for _, triangle := range t.tris { + // TODO: Find a way to quickly skip this triangle (or a way to iterate a subset of triangles) + testPointToTriangle := p.Sub(triangle.V[0]) + triNormal := triangle.Normal() + signedDistanceToTriPlane := triNormal.Dot(testPointToTriangle) + // Take this triangle as the source of truth if the projection of the point on the triangle is the closest + distToTri, _ := stlPointToTriangleDistSq(p, triangle) + if distToTri < closestTriangle { + closestTriangle = distToTri + signedDistanceResult = signedDistanceToTriPlane + } } - //log.Println("---") - ret := sigmoidScaled(retAccum / retWeight) - //log.Println(ret) - return ret -} - -func sigmoidScaled(x float64) float64 { - return 2/(1+math.Exp(-x)) - 1 // [-1, 1] + return signedDistanceResult } -func (m *triModelSdf) BoundingBox() sdf.Box3 { - return m.bb +func (t *triMeshSdf) BoundingBox() sdf.Box3 { + return t.bb } // ImportTriMesh converts a triangle-based mesh into a SDF3 surface. -func ImportTriMesh(tris chan *render.Triangle3, resolution float64, numSamples int) sdf.SDF3 { - m := &triModelSdf{ - samplesToTriangles: kdtree.New(nil), +// +// It is recommended to cache its values at setup time by using sdf.VoxelSdf. +// +// WARNING: It will only work on non-intersecting closed-surface(s) meshes. +// NOTE: Fix using blender for intersecting surfaces: Edit mode > P > By loose parts > Add boolean modifier to join them +func ImportTriMesh(tris []*render.Triangle3) sdf.SDF3 { + m := &triMeshSdf{ + tris: tris, bb: sdf.Box3{ Min: sdf.V3{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64}, Max: sdf.V3{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64}, }, - numSamples: numSamples, } - for triangle := range tris { - // Register the triangle's centroid associated to this triangle for later evaluation - triCentroid := triangle.V[0].Add(triangle.V[1]).Add(triangle.V[2]).DivScalar(3) - m.samplesToTriangles.Insert(points.NewPoint([]float64{triCentroid.X, triCentroid.Y, triCentroid.Z}, triangle)) - // Update the bounding box + + // Compute the bounding box + for _, triangle := range tris { for _, vertex := range triangle.V { - //vertex2 := vertex.MulScalar(0.99).Add(triCentroid.MulScalar(0.01)) - //m.samplesToTriangles.Insert(points.NewPoint([]float64{vertex2.X, vertex2.Y, vertex2.Z}, triangle)) m.bb = m.bb.Include(vertex) } - //// Sample inside the triangle (detail parameter) to avoid problems, - //// e.g. large triangles that are closer than small triangles that have closer vertices - //edge1 := triangle.V[1].Sub(triangle.V[0]) - //edge2 := triangle.V[2].Sub(triangle.V[0]) - //edgeLen := sdf.V2{X: edge1.Length(), Y: edge2.Length()} - //eps := 1e-3 - //uv := sdf.V2{} - //for uv.X = eps; uv.X < edgeLen.X-eps; uv.X += resolution { - // for uv.Y = eps; uv.Y < edgeLen.Y-eps; uv.Y += resolution { - // // Create the "uniform" sample point - // samplePoint := triangle.V[0].Add(edge1.MulScalar(uv.X)).Add(edge2.MulScalar(uv.Y)) - // // Ignore this sample if it lies outside the triangle - // if !stlPointInTriangle(samplePoint, triangle.V[0], triangle.V[1], triangle.V[2], 0) { // TODO: Efficiency? - // continue - // } - // // Add the "uniform" triangle sample to the tree - // m.samplesToTriangles.Insert(points.NewPoint([]float64{samplePoint.X, samplePoint.Y, samplePoint.Z}, triangle)) - // } - //} } if !m.bb.Contains(m.bb.Min) { // Return a valid bounding box if no vertices are found in the mesh m.bb = sdf.Box3{} // Empty box centered at {0,0,0} } - m.bb = m.bb.ScaleAboutCenter(1 + 1e-12) // Avoids missing faces due to inaccurate math operations. + //m.bb = m.bb.ScaleAboutCenter(1 + 1e-12) // Avoids missing faces due to inaccurate math operations. + return m } -func stlSameSide(p1, p2, a, b sdf.V3, tol float64) bool { - cp1 := b.Sub(a).Cross(p1.Sub(a)) - cp2 := b.Sub(a).Cross(p2.Sub(a)) - return cp1.Dot(cp2) >= -tol +//----------------------------------------------------------------------------- + +func stlPointToTriangleDistSq(p sdf.V3, triangle *render.Triangle3) (float64, bool /* falls outside? */) { + // Compute the closest point + closest, fallsOutside := stlClosestTrianglePointTo(p, triangle) + // Compute distance to the closest point + closestToP := p.Sub(closest) + distance := closestToP.Length2() + // Solve influence (distance) ties, by prioritizing triangles with normals more aligned to `closestToP`. + // This should fix ghost triangle extensions and smooth the field over sharp angles. + if fallsOutside { // <-- This is an optimization, as others have 0 extra influence in this step + distance *= 1 + (1-math.Abs(closestToP.Normalize().Dot(triangle.Normal())))*stlEpsilon + } + //log.Println(distance, closestToP.Normalize().Dot(triangle.Normal())) + return distance, fallsOutside +} + +// https://stackoverflow.com/a/47505833 +func stlClosestTrianglePointTo(p sdf.V3, triangle *render.Triangle3) (sdf.V3, bool /* falls outside? */) { + edgeAbDelta := triangle.V[1].Sub(triangle.V[0]) + edgeCaDelta := triangle.V[0].Sub(triangle.V[2]) + edgeBcDelta := triangle.V[2].Sub(triangle.V[1]) + + // The closest point may be a vertex + uab := stlEdgeProject(triangle.V[0], edgeAbDelta, p) + uca := stlEdgeProject(triangle.V[2], edgeCaDelta, p) + if uca > 1 && uab < 0 { + return triangle.V[0], true + } + ubc := stlEdgeProject(triangle.V[1], edgeBcDelta, p) + if uab > 1 && ubc < 0 { + return triangle.V[1], true + } + if ubc > 1 && uca < 0 { + return triangle.V[2], true + } + + // The closest point may be on an edge + triNormal := triangle.Normal() + planeAbNormal := triNormal.Cross(edgeAbDelta) + planeBcNormal := triNormal.Cross(edgeBcDelta) + planeCaNormal := triNormal.Cross(edgeCaDelta) + if uab >= 0 && uab <= 1 && !stlPlaneIsAbove(triangle.V[0], planeAbNormal, p) { + return stlEdgePointAt(triangle.V[0], edgeAbDelta, uab), true + } + if ubc >= 0 && ubc <= 1 && !stlPlaneIsAbove(triangle.V[1], planeBcNormal, p) { + return stlEdgePointAt(triangle.V[1], edgeBcDelta, ubc), true + } + if uca >= 0 && uca <= 1 && !stlPlaneIsAbove(triangle.V[2], planeCaNormal, p) { + return stlEdgePointAt(triangle.V[2], edgeCaDelta, uca), true + } + + // The closest point is in the triangle so project to the plane to find it + return stlPlaneProject(triangle.V[0], triNormal, p), false +} + +func stlEdgeProject(edge1, edgeDelta, p sdf.V3) float64 { + return p.Sub(edge1).Dot(edgeDelta) / edgeDelta.Length2() } -func stlPointInTriangle(p, a, b, c sdf.V3, tol float64) bool { - return stlSameSide(p, a, b, c, tol) && stlSameSide(p, b, a, c, tol) && stlSameSide(p, c, a, b, tol) +func stlEdgePointAt(edge1, edgeDelta sdf.V3, t float64) sdf.V3 { + return edge1.Add(edgeDelta.MulScalar(t)) } -// ImportSTL converts an STL model into a SDF3 surface. -func ImportSTL(reader io.ReadSeeker, resolution float64, numSamples int) (sdf.SDF3, error) { +func stlPlaneIsAbove(anyPoint, normal, testPoint sdf.V3) bool { + return normal.Dot(testPoint.Sub(anyPoint)) > 0 +} + +func stlPlaneProject(anyPoint, normal, testPoint sdf.V3) sdf.V3 { + v := testPoint.Sub(anyPoint) + d := normal.Dot(v) + p := testPoint.Sub(normal.MulScalar(d)) + return p +} + +//----------------------------------------------------------------------------- + +// ImportSTL converts an STL model into a SDF3 surface. See ImportTriMesh. +func ImportSTL(reader io.ReadSeeker) (sdf.SDF3, error) { mesh, err := stl.ReadAll(reader) if err != nil { return nil, err } - triChan := make(chan *render.Triangle3, 64) // Buffer some triangles and send in batches if scheduler prefers it - go func() { - for _, triangle := range mesh.Triangles { - tri := &render.Triangle3{} - for i, vertex := range triangle.Vertices { - tri.V[i] = stlToV3(vertex) - } - triChan <- tri + tris := make([]*render.Triangle3, 0) // Buffer some triangles and send in batches if scheduler prefers it + for _, triangle := range mesh.Triangles { + tri := &render.Triangle3{} + for i, vertex := range triangle.Vertices { + tri.V[i] = sdf.V3{X: float64(vertex[0]), Y: float64(vertex[1]), Z: float64(vertex[2])} } - close(triChan) - }() - return ImportTriMesh(triChan, resolution, numSamples), nil -} - -func stlToV3(v stl.Vec3) sdf.V3 { - return sdf.V3{X: float64(v[0]), Y: float64(v[1]), Z: float64(v[2])} + tris = append(tris, tri) + } + return ImportTriMesh(tris), nil } diff --git a/sdf/voxel.go b/sdf/voxel.go new file mode 100644 index 000000000..0c20659a8 --- /dev/null +++ b/sdf/voxel.go @@ -0,0 +1,89 @@ +//----------------------------------------------------------------------------- +/* + +Voxel-based cache to remove deep SDF3 hierarchies at setup and speed up evaluation + +*/ +//----------------------------------------------------------------------------- + +package sdf + +// VoxelSdf is the SDF that represents a pre-computed voxel-based SDF3. +//It can be used as a cache, or for smoothing. +// +// CACHE: +// It can be used to speed up all evaluations required by the surface mesher at the cost of scene setup time and accuracy. +// +// SMOOTHING (meshCells <<< renderer's meshCells): +// It performs trilinear mapping for inner values and may be used as a cache for any other SDF, losing some accuracy. +// +// WARNING: It may lose sharp features, even if meshCells is high. +type VoxelSdf struct { + // voxelCorners are the values of this SDF in each voxel corner + voxelCorners map[V3i]float64 // TODO: Octree + k-d tree to simplify/reduce memory consumption + speed-up access? + // bb is the bounding box. + bb Box3 + // Number of voxelCorners to consider + numVoxels V3i +} + +// NewVoxelSDF see VoxelSdf. This populates the whole cache from the given SDF. The progress listener may be nil. +func NewVoxelSDF(s SDF3, meshCells int, progress chan float64) SDF3 { + bb := s.BoundingBox() // TODO: Use default code to avoid duplication + bbSize := bb.Size() + resolution := bbSize.MaxComponent() / float64(meshCells) + cells := bbSize.DivScalar(resolution).ToV3i() + + voxelCorners := map[V3i]float64{} + voxelCornerIndex := V3i{} + for voxelCornerIndex[0] = 0; voxelCornerIndex[0] <= cells[0]; voxelCornerIndex[0]++ { + for voxelCornerIndex[1] = 0; voxelCornerIndex[1] <= cells[1]; voxelCornerIndex[1]++ { + for voxelCornerIndex[2] = 0; voxelCornerIndex[2] <= cells[2]; voxelCornerIndex[2]++ { + voxelCorner := bb.Min.Add(bbSize.Mul(voxelCornerIndex.ToV3()).Div(cells.ToV3())) + voxelCorners[voxelCornerIndex] = s.Evaluate(voxelCorner) + } + } + if progress != nil { + progress <- float64(voxelCornerIndex[0]) / float64(cells[0]) + } + } + + return &VoxelSdf{ + voxelCorners: voxelCorners, + bb: bb, + numVoxels: cells, + } +} + +func (m *VoxelSdf) Evaluate(p V3) float64 { + // Find the voxel's {0,0,0} corner quickly and compute p's displacement + voxelSize := m.bb.Size().Div(m.numVoxels.ToV3()) + voxelStartIndex := p.Sub(m.bb.Min).Div(voxelSize).ToV3i() + voxelStart := m.bb.Min.Add(voxelSize.Mul(voxelStartIndex.ToV3())) + d := p.Sub(voxelStart).Div(voxelSize) // [0, 1) for each dimension + // Get the values at the voxel's corners + c000 := m.voxelCorners[voxelStartIndex] + c001 := m.voxelCorners[voxelStartIndex.Add(V3i{0, 0, 1})] + c010 := m.voxelCorners[voxelStartIndex.Add(V3i{0, 1, 0})] + c011 := m.voxelCorners[voxelStartIndex.Add(V3i{0, 1, 1})] + c100 := m.voxelCorners[voxelStartIndex.Add(V3i{1, 0, 0})] + c101 := m.voxelCorners[voxelStartIndex.Add(V3i{1, 0, 1})] + c110 := m.voxelCorners[voxelStartIndex.Add(V3i{1, 1, 0})] + c111 := m.voxelCorners[voxelStartIndex.Add(V3i{1, 1, 1})] + // Perform trilinear interpolation over the voxel's corners + // - 4 linear interpolations + c00 := c000*(1-d.X) + c100*d.X + c01 := c001*(1-d.X) + c101*d.X + c10 := c010*(1-d.X) + c110*d.X + c11 := c011*(1-d.X) + c111*d.X + // - 2 bilinear interpolations + c0 := c00*(1-d.Y) + c10*d.Y + c1 := c01*(1-d.Y) + c11*d.Y + // - 1 trilinear interpolation + c := c0*(1-d.Z) + c1*d.Z + return c +} + +func (m *VoxelSdf) BoundingBox() Box3 { + return m.bb +} From 45eb220717438a224f7ad25693591121ab79fe05 Mon Sep 17 00:00:00 2001 From: Yeicor Date: Sat, 1 Jan 2022 10:46:02 +0100 Subject: [PATCH 3/4] Update SHA1 for new example --- examples/monkey_hat/SHA1SUM | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/monkey_hat/SHA1SUM b/examples/monkey_hat/SHA1SUM index 688b1056a..5c0a1b17a 100644 --- a/examples/monkey_hat/SHA1SUM +++ b/examples/monkey_hat/SHA1SUM @@ -1,2 +1 @@ -6712d16e1047f5281e2968ee4543a9aaac46fb40 pool2.stl -e2d04358f90044d6b5f73cd17df3e9c55e70267c pool1.stl +d717263395e88df0a1f878b6e51021b02f0656ba monkey-out.stl From ccd3337b79d871c967e7145b6396b96ffa3923ea Mon Sep 17 00:00:00 2001 From: Yeicor Date: Sat, 8 Jan 2022 19:38:10 +0100 Subject: [PATCH 4/4] Heavily optimized STL importer --- examples/monkey_hat/main.go | 5 +-- go.mod | 1 + go.sum | 2 + obj/stl.go | 80 ++++++++++++++++++++++++++----------- sdf/voxel.go | 16 ++++---- 5 files changed, 70 insertions(+), 34 deletions(-) diff --git a/examples/monkey_hat/main.go b/examples/monkey_hat/main.go index 590b9f6d1..0b23d2be4 100644 --- a/examples/monkey_hat/main.go +++ b/examples/monkey_hat/main.go @@ -28,11 +28,10 @@ func monkeyWithHat() sdf.SDF3 { } } // - Create the SDF from the mesh (a modified Suzanne from Blender with 366 faces) - monkeyImported, err := obj.ImportSTL(file) + monkeyImported, err := obj.ImportSTL(file, 20, 3, 5) if err != nil { panic(err) } - //monkeyImported = sdf.NewVoxelSDF(monkeyImported, 64) // HAT hatHeight := 0.5 @@ -55,7 +54,7 @@ func monkeyWithHat() sdf.SDF3 { // It also smooths the mesh a little using trilinear interpolation. // It is actually slower for this mesh (unless meshCells <<< renderer's meshCells), but should be faster for // more complex meshes (with more triangles) or SDF3 hierarchies that take longer to evaluate. - monkeyHat = sdf.NewVoxelSDF(monkeyHat, 64, nil) // Use 32 for harder smoothing demo + monkeyHat = sdf.NewVoxelSDF3(monkeyHat, 64, nil) // Use 32 for harder smoothing demo return monkeyHat } diff --git a/go.mod b/go.mod index 4520fdd0c..95cebcb8e 100644 --- a/go.mod +++ b/go.mod @@ -4,6 +4,7 @@ go 1.13 require ( github.com/ajstarks/svgo v0.0.0-20200725142600-7a3c8b57fecb + github.com/dhconnelly/rtreego v1.1.0 github.com/golang/freetype v0.0.0-20170609003504-e2365dfdc4a0 github.com/hschendel/stl v1.0.4 github.com/llgcode/draw2d v0.0.0-20200930101115-bfaf5d914d1e diff --git a/go.sum b/go.sum index 2b9d1924e..41bf66bb8 100644 --- a/go.sum +++ b/go.sum @@ -7,6 +7,8 @@ github.com/ajstarks/svgo v0.0.0-20200725142600-7a3c8b57fecb/go.mod h1:K08gAheRH3 github.com/boombuler/barcode v1.0.0/go.mod h1:paBWMcWSl3LHKBqUq+rly7CNSldXjb2rDl3JlRe0mD8= github.com/davecgh/go-spew v1.1.0 h1:ZDRjVQ15GmhC3fiQ8ni8+OwkZQO4DARzQgrnXU1Liz8= github.com/davecgh/go-spew v1.1.0/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38= +github.com/dhconnelly/rtreego v1.1.0 h1:ejMaqN03N1s6Bdg6peGkNgBnYYSBHzcK8yhSPCB+rHE= +github.com/dhconnelly/rtreego v1.1.0/go.mod h1:SDozu0Fjy17XH1svEXJgdYq8Tah6Zjfa/4Q33Z80+KM= github.com/fogleman/gg v1.2.1-0.20190220221249-0403632d5b90/go.mod h1:R/bRT+9gY/C5z7JzPU0zXsXHKM4/ayA+zqcVNZzPa1k= github.com/fogleman/gg v1.3.0/go.mod h1:R/bRT+9gY/C5z7JzPU0zXsXHKM4/ayA+zqcVNZzPa1k= github.com/go-fonts/dejavu v0.1.0/go.mod h1:4Wt4I4OU2Nq9asgDCteaAaWZOV24E+0/Pwo0gppep4g= diff --git a/obj/stl.go b/obj/stl.go index e808bb9db..3cfde91da 100644 --- a/obj/stl.go +++ b/obj/stl.go @@ -11,6 +11,7 @@ package obj import ( "github.com/deadsy/sdfx/render" "github.com/deadsy/sdfx/sdf" + "github.com/dhconnelly/rtreego" "github.com/hschendel/stl" "io" "math" @@ -19,22 +20,21 @@ import ( //----------------------------------------------------------------------------- type triMeshSdf struct { - tris []*render.Triangle3 - bb sdf.Box3 + rtree *rtreego.Rtree + numNeighbors int + bb sdf.Box3 } -const stlEpsilon = 1e-12 +const stlEpsilon = 1e-1 func (t *triMeshSdf) Evaluate(p sdf.V3) float64 { - if !t.bb.Contains(p) { // Fast exit - // Length to surface is at least distance to bounding box - return t.bb.Include(p).Size().Sub(t.bb.Size()).Length() - } // Check all triangle distances signedDistanceResult := 1. closestTriangle := math.MaxFloat64 - for _, triangle := range t.tris { - // TODO: Find a way to quickly skip this triangle (or a way to iterate a subset of triangles) + // Quickly skip checking most triangles by only checking the N closest neighbours (AABB based) + neighbors := t.rtree.NearestNeighbors(t.numNeighbors, stlToPoint(p)) + for _, neighbor := range neighbors { + triangle := neighbor.(*stlTriangle).Triangle3 testPointToTriangle := p.Sub(triangle.V[0]) triNormal := triangle.Normal() signedDistanceToTriPlane := triNormal.Dot(testPointToTriangle) @@ -52,15 +52,22 @@ func (t *triMeshSdf) BoundingBox() sdf.Box3 { return t.bb } -// ImportTriMesh converts a triangle-based mesh into a SDF3 surface. +// ImportTriMesh converts a triangle-based mesh into a SDF3 surface. minChildren and maxChildren are parameters that can +// affect the performance of the internal data structure (3 and 5 are a good default; maxChildren >= minChildren > 0). +// +// WARNING: Setting a low numNeighbors will consider many fewer triangles for each evaluated point, greatly speeding up +// the algorithm. However, if the count of triangles is too low artifacts will appear on the surface (triangle +// continuations). Setting this value to MaxInt is extremely slow but will provide correct results, so choose a value +// that works for your model. // -// It is recommended to cache its values at setup time by using sdf.VoxelSdf. +// It is recommended to cache (and/or smooth) its values by using sdf.VoxelSdf3. // // WARNING: It will only work on non-intersecting closed-surface(s) meshes. // NOTE: Fix using blender for intersecting surfaces: Edit mode > P > By loose parts > Add boolean modifier to join them -func ImportTriMesh(tris []*render.Triangle3) sdf.SDF3 { +func ImportTriMesh(tris chan *render.Triangle3, numNeighbors, minChildren, maxChildren int) sdf.SDF3 { m := &triMeshSdf{ - tris: tris, + rtree: nil, + numNeighbors: numNeighbors, bb: sdf.Box3{ Min: sdf.V3{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64}, Max: sdf.V3{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64}, @@ -68,7 +75,9 @@ func ImportTriMesh(tris []*render.Triangle3) sdf.SDF3 { } // Compute the bounding box - for _, triangle := range tris { + bulkLoad := make([]rtreego.Spatial, 0) + for triangle := range tris { + bulkLoad = append(bulkLoad, &stlTriangle{Triangle3: triangle}) for _, vertex := range triangle.V { m.bb = m.bb.Include(vertex) } @@ -77,6 +86,7 @@ func ImportTriMesh(tris []*render.Triangle3) sdf.SDF3 { m.bb = sdf.Box3{} // Empty box centered at {0,0,0} } //m.bb = m.bb.ScaleAboutCenter(1 + 1e-12) // Avoids missing faces due to inaccurate math operations. + m.rtree = rtreego.NewTree(3, minChildren, maxChildren, bulkLoad...) return m } @@ -158,19 +168,43 @@ func stlPlaneProject(anyPoint, normal, testPoint sdf.V3) sdf.V3 { //----------------------------------------------------------------------------- +type stlTriangle struct { + *render.Triangle3 +} + +func (s *stlTriangle) Bounds() *rtreego.Rect { + bounds := sdf.Box3{Min: s.V[0], Max: s.V[0]} + bounds = bounds.Include(s.V[1]) + bounds = bounds.Include(s.V[2]) + points, err := rtreego.NewRectFromPoints(stlToPoint(bounds.Min), stlToPoint(bounds.Max)) + if err != nil { + panic(err) // Implementation error + } + return points +} + +func stlToPoint(v3 sdf.V3) rtreego.Point { + return rtreego.Point{v3.X, v3.Y, v3.Z} +} + +//----------------------------------------------------------------------------- + // ImportSTL converts an STL model into a SDF3 surface. See ImportTriMesh. -func ImportSTL(reader io.ReadSeeker) (sdf.SDF3, error) { +func ImportSTL(reader io.ReadSeeker, numNeighbors, minChildren, maxChildren int) (sdf.SDF3, error) { mesh, err := stl.ReadAll(reader) if err != nil { return nil, err } - tris := make([]*render.Triangle3, 0) // Buffer some triangles and send in batches if scheduler prefers it - for _, triangle := range mesh.Triangles { - tri := &render.Triangle3{} - for i, vertex := range triangle.Vertices { - tri.V[i] = sdf.V3{X: float64(vertex[0]), Y: float64(vertex[1]), Z: float64(vertex[2])} + tris := make(chan *render.Triangle3, 128) // Buffer some triangles and send in batches if scheduler prefers it + go func() { + for _, triangle := range mesh.Triangles { + tri := &render.Triangle3{} + for i, vertex := range triangle.Vertices { + tri.V[i] = sdf.V3{X: float64(vertex[0]), Y: float64(vertex[1]), Z: float64(vertex[2])} + } + tris <- tri } - tris = append(tris, tri) - } - return ImportTriMesh(tris), nil + close(tris) + }() + return ImportTriMesh(tris, numNeighbors, minChildren, maxChildren), nil } diff --git a/sdf/voxel.go b/sdf/voxel.go index 0c20659a8..2601ccd53 100644 --- a/sdf/voxel.go +++ b/sdf/voxel.go @@ -1,14 +1,14 @@ //----------------------------------------------------------------------------- /* -Voxel-based cache to remove deep SDF3 hierarchies at setup and speed up evaluation +Voxel-based cache/smoothing to remove deep SDF2/SDF3 hierarchies and speed up evaluation */ //----------------------------------------------------------------------------- package sdf -// VoxelSdf is the SDF that represents a pre-computed voxel-based SDF3. +// VoxelSdf3 is the SDF that represents a pre-computed voxel-based SDF3. //It can be used as a cache, or for smoothing. // // CACHE: @@ -18,7 +18,7 @@ package sdf // It performs trilinear mapping for inner values and may be used as a cache for any other SDF, losing some accuracy. // // WARNING: It may lose sharp features, even if meshCells is high. -type VoxelSdf struct { +type VoxelSdf3 struct { // voxelCorners are the values of this SDF in each voxel corner voxelCorners map[V3i]float64 // TODO: Octree + k-d tree to simplify/reduce memory consumption + speed-up access? // bb is the bounding box. @@ -27,8 +27,8 @@ type VoxelSdf struct { numVoxels V3i } -// NewVoxelSDF see VoxelSdf. This populates the whole cache from the given SDF. The progress listener may be nil. -func NewVoxelSDF(s SDF3, meshCells int, progress chan float64) SDF3 { +// NewVoxelSDF3 see VoxelSdf3. This populates the whole cache from the given SDF. The progress listener may be nil. +func NewVoxelSDF3(s SDF3, meshCells int, progress chan float64) SDF3 { bb := s.BoundingBox() // TODO: Use default code to avoid duplication bbSize := bb.Size() resolution := bbSize.MaxComponent() / float64(meshCells) @@ -48,14 +48,14 @@ func NewVoxelSDF(s SDF3, meshCells int, progress chan float64) SDF3 { } } - return &VoxelSdf{ + return &VoxelSdf3{ voxelCorners: voxelCorners, bb: bb, numVoxels: cells, } } -func (m *VoxelSdf) Evaluate(p V3) float64 { +func (m *VoxelSdf3) Evaluate(p V3) float64 { // Find the voxel's {0,0,0} corner quickly and compute p's displacement voxelSize := m.bb.Size().Div(m.numVoxels.ToV3()) voxelStartIndex := p.Sub(m.bb.Min).Div(voxelSize).ToV3i() @@ -84,6 +84,6 @@ func (m *VoxelSdf) Evaluate(p V3) float64 { return c } -func (m *VoxelSdf) BoundingBox() Box3 { +func (m *VoxelSdf3) BoundingBox() Box3 { return m.bb }