From 07435b1378e94983965977483974740268fcd391 Mon Sep 17 00:00:00 2001 From: Kuldeep Kumar Date: Mon, 5 Dec 2016 16:53:34 +0100 Subject: [PATCH] Varifold with EP and MSM dec 5 2016 Commit --- CMakeLists.txt | 29 ++ MSmeasure_GFA_Sub1_test50fib | Bin 0 -> 3000 bytes Tracts_Sub1_VTK_test50fib.fib | Bin 0 -> 14590 bytes compute_gramiam_varifolds.cpp | 375 ++++++++++++++++++++ compute_varifold_using_EP_and_MSmeasure.cpp | 80 +++-- 5 files changed, 454 insertions(+), 30 deletions(-) create mode 100644 CMakeLists.txt create mode 100644 MSmeasure_GFA_Sub1_test50fib create mode 100644 Tracts_Sub1_VTK_test50fib.fib create mode 100644 compute_gramiam_varifolds.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..ffe5f08 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,29 @@ + +cmake_minimum_required(VERSION 2.8) + +PROJECT(compute_varifold_using_EP_and_MSmeasure) + +find_package(VTK REQUIRED) +include(${VTK_USE_FILE}) + +find_package(OpenMP) +if (OPENMP_FOUND) + set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") +endif() + +# Find ITK. +FIND_PACKAGE(ITK REQUIRED) +IF(ITK_FOUND) + INCLUDE(${ITK_USE_FILE}) +ENDIF(ITK_FOUND) + +INCLUDE_DIRECTORIES (/usr/include/eigen3) + +add_executable(compute_varifold_using_EP_and_MSmeasure compute_varifold_using_EP_and_MSmeasure.cpp) + +if(VTK_LIBRARIES) + target_link_libraries(compute_varifold_using_EP_and_MSmeasure ${VTK_LIBRARIES} ITKCommon ) +else() + target_link_libraries(compute_varifold_using_EP_and_MSmeasure vtkHybrid vtkWidgets) +endif() diff --git a/MSmeasure_GFA_Sub1_test50fib b/MSmeasure_GFA_Sub1_test50fib new file mode 100644 index 0000000000000000000000000000000000000000..c7a32be34852074cb8b8a47672a0b0d39e5dc75d GIT binary patch literal 3000 zcmWNS2|E=C1BLBL=yH{%dn2I~AzP($rp=a`MpLP#GA5-Wno5$UMYJfD15Q`~%lup6?&FW($2~Y@8i;U>3Ur+Ngo_Kk;n@_A;Kkk$yV{_B z@G00_K8dbIH+=9=V}xEC6ePwF9Iu>e?Gu3G;Xb%K_9KFn^he^i66E67(a>q6_xTG0+Gf#7{>4z*;Q^t3=~csVpCE zvEjt-03-&h_G>)U`+ zd+x*9w-9Do{g@DB!|cJKxLL3qc0+RDVG@ntZ8vdbM=C<}2l6D}!`Kz>XgIST+@Zw} zAyVvck*C=bRVG&0ak}qoXeZfI%d!pEzA5wEzGirNOR)XNZE&p;TisQdaJUP$x@XZ4 zt4og(Kgd{{^RR6a79|hjy({zCoD+$*V+NGko`PtNFXC;BaV15bHxkNWG4dx)EfZ4W z_6i#urlR-mC8V$J#XqZBk@w~WVs8lP9tPfPHsGy-2xyA8;h*#AupEB_9o41q$;-g7 zt#6@qP?c52;V`^@84vFLK)|5}Y$g2H(eZ&0ZYtnZSTQvf$OB$KxB;F)>X{w=u^d zsZB>)vI9Q<^BxE03i;eS3wu4ZcraunV;5h9wOkRh>y_!9uf{2&DvYiAg5H=ZEI9WN zcC&urmUtuH5gE{V_$2zB(WLJdDJF#c#0U*-9^Lr~E646f;4cbn*LZWjL3Y-+_VZRJ?tH=eI#MJ&IZhDnX{sdrMHSIAyk$@JC4op6igE+}j_V2%S z^1h`h+(RIkoy%~IArH=b{d?8fy#82c2^-{v5?V+ z_aeGC1*&;!43rMxgg_nI-W$rh54NJf#hJk!)rhU1#n8kQ%sJx8VL?g^oZb&VM{P!k zpQ20LgNC66kjw_Vi;iGgh%T{NpU;&?aQL@%XzZ%N*RW66bZ`NcURR=Sqy)JxzcE~K z6PK4p;;UH|?(FNv8f0L~y&6Ow_U77gX;^qx$a8sF_%g_bHFb(yx=)#F=5|5dz>Bj| zE!eZL5GyQnSa`Y~g=@5U^Zg7ATa$%|P7}KH38@%kNb~X5EHMm#mCI5t-P#G=(Z+n% z^co-kTELX9TI89ivC>|P20InFw6Y1KGpbRkVnM+LDf~|f==noP{{khN-P?@-rE+|` zHi&=BoXmjJAO_zm$hhf@4}C}WQ9Mb8x(s*kc{Mzk(prs`Av zh%pm)ti_#|^EviERXUc*Q8(LwaUL$b-|!jfrpkOAS%D>jYRG-dgRN;F>bHnESs=r{ zR9O}#8gOHr97kV@#f%}jXv!|bzw>&KGG}0S&ezcyWX;Se+gR!uiW6OB(AcX??a{qB zscBCQXDg0;q|AtNJ+Aku#s8WvVaF|9c1$s(f1VE2ZW;0869?|jmEo%eEs&TB>0Hx; zpHWwl6FHe*^G*4^L5~4uVm^O9f?@xL!flr!qtxFc*Y7PlO>{WslqG94deO3LC^x=q z!Q2yPFh=h>=Elk}k3xF8cH>`NW%^qWW5}m$gv-WZM8ZQP^~=)vc^y*i2YB>ViX0y! zCf#*qgz0>!IJM)=TQUE*p~wbD8w%76>7*ym_*g^kDXznV^|^?`2zt+(%)y;XocBVF zOSkJYH(?ZWb!zc!u@bMge#EHvirk~?%-Ty9Ow-Zib#-0pbWY@!ZGlM8AI5EUC8*GN zgZoB(n6|JTG}dLBuNH4v>F`}vIOc@i#M}W#iN~t4{6`b|I*d5gtsjR5+_>=2BeZO& zLf{|c_;Gy)W=0utMt31fGlZPspusm^R-q^*8}l6oahzvAJnoNTX_7wupOnF7qX9=P ztHkhDF}ueZ@p-g8BLbCJ6{f%!D|!(cHif>$Vz#YOU{78R$~H{li6fQ_Qyk6YM0u7p zXtH~uEPeZTbGQ}LrKUAO$pkh672bIPQfk>%D8Ay^k@%Mf-Xb4JO}MXA93wB z2_|NnQq!p)l^IV^(N>F`tqRQEC!(2CCpJE9hNwk?WBCHg6q~a{7>0#0X^2#Bg`9dm zX1c30d7hXj1&6c(n(e4y{o7IC2YR!lr@LtAOH_ChUprXAN zD{qLnGf|JaHU@Ou{S#JJ$~?cwfm@S9VHz@+p3*y5a%~R(?oP#jtG{D`_D}}#Goq5D zh+Y+i)7K6lQLzw#`NdGy&&A)l#*Ewi8YQ)_F=yIKJdG9dMXQLOd5yS|aRG_-7WD2? z|Dm?n}1C}h3Vf;Nq9-0?|rE)4%pB0J* zsVeo}?7{FNMLN7~z%0ReSZ!8dS>_@58T;U0_ghHSzKZJRc1(^6=luY$|50vxP= zjB`ORacf-z?CvZ;w`LB$Bu_=DP=$B0cENsc5p3NGP$TC7uN5T-$=DAYS2u{JreU{C z+(3668g=$AL?|RPm1K$xnWscThBD8J%t_`5oqZ{iDHMei z6_I(!SU%f&`{Da9eAjxldaV1oRy`Z8b)Wm}>psWEa#q|V%YcZuh#nJS##)Arojzki z>=aAaj?Nl?zYse!ZtBc94eZqq1`O2*fKV9NJ=8KJxX;kPStBI4SHIAKmR;SQEu&*% zBjPl`&T%RxkeXk1lmY|~cg2{EcO=E>-jGmQfT3rO(FYxyWwY4 z7@F>UPSIwWsBq>0njCmZsgZ|JW$-REIuk-!?t!SCl7RX@&(RGWh6Z~ZQLA+Xy&CWa z4HiY9%7W#hnfF!Hb^M45yABBJ$RgDB>;T%E9E8%VjlwbVru58t2EBclD4g%ONfQqK zrpni?gsYz8<0g}Fh2&|j!W?c_beroK=ppEM@dly{~xq5`B=n*x84v*X`&rOg8ICh}1``t(tP5~$$kl5%z2DXeGJcE`^s0RS7AhNUAXtql@#lYaky7IxHhdHS>H3j&?DF2^5^g5 zJ$5|ywf+XDg6B~Gwze3czY>mCG*DRW8T5&K0y!`2XkaH@>{1a6H@6s4&$S)Vbx#vK zi*8Pxn=H}I{1JSOoByB(`QGm;k+`=tEzceRxmvGLDHF=F)%O4l7-5C>O69Pp%XNsn;DVjY<6v`^7A$V8MZf%M zkfw15R!q2xeFr14iES@Ffts~7@IAx{t&;tr z{9RA@I9Nt=T`#z8U;*`>253068ie~Ys0nmHolFfleJ2T?j@Cl8A9ZlB#Z9PnbwNd~ zG&q_X3m@kXhhKB{!`Yl>&@?p?e!o+I{J<1c{+))3B_rVO&RM7)*aKB6pF(8|f7CLm zM6Ilq@aE}w)KNK#EpFd~p9N8Na&DvcC@-?QycL5V`N|*ru-ujbKckn#+}tFU+O4NsW62t$6CTk>@Ej-Rl`)j8YsA+Ae*^JQdN2y)ZYo9#h%VWxz!Pj zv~i;ug@Z+NX){J!jHBq@U4_c#Je=5a2MxA)CDhE;;N*jwsplbIp?;|+PSbr$o-b>K z=G;CQ_oM@LRB0o$L;Y~hFJrRZ-b3iT-GK|&XHy%$i$d>a7B2C6Os&uV6fFjI$ECLl z$i}5ow5Z;L%VSN*?&v(x;?h}MX=Y1pH0^{zvrU+M#E^^wP#8FS;acm?q;|kV=;^J- z_3w>fTdOXVt>cU;4o~6W;oWp{*=N*iJ_fQa*V2hDQ_w=?1f2BfPscS&(c#E1I97I@ zj%Kw+cZ27!?~*wk3baG%;8)m^P(vA2J{WW{8dh~LrL==DFl@RJ%+In+7-trI0;@HDXlGG1F22(Nv{HM~4$XI%c(gfawi`)nn>J#y zb3NFctRgHrgzI#jM6cIRq}z_plpMQS`1p;#d~}uzt-gFsIO%Mc?-hQa710}nlS>Qg zcPXAcmMDmJwz;%TW2sz6Wx@jI(d{U0*^9mjqb42E%qY@r1{#+Kpj-Qv!q4b~aF@=YTkllS6O2Sxqs{2Lv4aRM zJ1#u7Y{kwxYs8S4AHp-K6FQp?7U2#3gy+UGbn4PiOz0LPJU1z!!^<==wLDUI_Hag< zMHj^EvueVld=a*~X(|>4y$~LXd1&^|Tr6#IO?0idMw3NR|pDfmwE)`u%OW>DRj@S~vTeuCkhEAs>T=v+yo9i|O7;#n! zr!IKx*MC6^NK9{zVaHeag|BdcZF_#Bug&_r@L)sO;X4Se#!S0>e^oBbIsXCGPi*u} zG>Zd^30L6j4KKeBuUbn!XU9YRjssFiA4!fsvl`w<50_Nkl*q=c0lw|sl$SBSfyUhY zhB~V5^5%vjTIL&%4*Fx{{O4)3aNkM{Fdi!(i!GCk!jNUhSf4{u2gJ2#0$ zrNvO$?mc~I8zK_M*FfpvH&i`AO~glig8Q?asI;pr7L|X3G9PU!TsmIFTdjwx`A#Id z4-rc$w!)jrD9S5}5lcTI{3t7I>Q`?Tkf6wyELb;(;g88Mk$)4Zf?_e8p+Fo;|}G>4E$m*~pX z60+#$0Q+`Y(yfl6WcV!*DjwXV+^l*sX|O~s^%y$Xq(kkmYN12uWlFt%kph%%pnqdL z%~z};w+agkANq?@io(gPOv2R#vCw5z2wV<&gabD{g`PInaISVDANyTlklP)|dfyBE zFRg=cqoHv4g(iBPnFrC@gCR3#1iG|$hw(OF;LO1$beP)&<8SQ-5iFx+!()gx{ss3| z)}h5r6Bs!`52}vbMC*1#q2D!v7iagOW9}dDHS>U1%k|K8+!ttH5D#ypmgr;o3rwR; z;Y0uZ7!Vr(&GwbTmnT24&+PZonh#a*v7;6aF`Xx^yQ3flt6I{YheQ2z0<@%1lNpqL zxK!?y^heU0Ur5_e7|PeYd3m3$RcX5OPCuI#jN3*WKB5jWbI{hM$w)e;<#j!)^ zMC1TEt@28~cz83NU0_I;_j&Uco49c1vPt~0>rRT1^H;JmF>bp7c zeV!Uv9llM)8P}jXQeV2F`IriC^?`dWUicMeDAAS9&Tw;aupH3#4P9xT0LAHD;XV^I0zT3ThfRhU)Bq2>o~ z+7z)|I5x~ceR~($w&bz!)viLrssdWG=9cJN*#S+;JJOP`2Sh*nF=+34ovhY}h`zIb zqqLmo=!=CyqbUlqcNYrJ^+$yE$|yM9MOnC1E*2VZM#6!0hmYFOpYPt!voztM?YOOFIK3(W7+zn5397JoCVZ!w4K&Yxu7S>U2!t(TfsQR-~ zv^zIn*v%RMPhwm}=ge-x-upb1KYk^=g^{ouUk_#Xt3|->8^U%)I+Q$lAo|=a7M3S1 zp`Tiq1zU#xWn3 zJ&eF+Gt{VekSBodXEeCoiA)vdNYS?0!fQurXJI@rH&N&>;Y)IlKZgpMn&$a_?sfD_hb(ZSXZd6S<-<5d=B)Izfgl@j7{@En98mfk2u)h zb^$zCm?al>*Q2@|1xn7DBdwjVm8x&Nrp=pPNV)sAQk8)drQb}D-u@g#70t>i^V)aN z``noxO?ye_Jj$W{p9xg9G>WK00QCHsK@T2wq@wy;5Y>M)m6j?~$&9WLzicu+_*qC5 zKVQMxCq7jEbs4?*Sqi)Q4yB6NH}t;x3miTEhMr~V(ziAr;KFr7s;$tVrY>4={fD5s zpc0`p1EKtkJv~#JDpaDk!%ufh%B)Mn$$HPFZO&=5Gba_tEZZlgE>fmls)$3{uaYF6 z!*uX{6ZUyBT?*jY^W?i*2o5!pb6zp!ozg@1?>nWQ7o6#a>ICeNlO^q(Sx7e>>dE67vPGAg_}gBn80tXX?E}&4?P;E!ooV^CnP_z78B8fS zM_wm~L(P=yknu#H^lh|Yb+_$s>jKF2Hru63KYqZwb8qF+3-z*G))AGr2Fmv{Zqs_F zL#R14pIU@j)2D57QCnk|{G&{w#vM7>@^_uQ_v0|yU1NvF35N3d1b_M0A}`chuP48o zG+J6P&K;VRgUHbQ3e2f2f{Lp#)Gh5I+@Imk&yidT9PNxw814`PKOE2|L zV6PXGXx8|x^r9&R2WniTQSWl-^^+tVVV6iDU)obc-aCvsJCFjcvV~Hmz?fx^$n|!L zP+R1RabYoJ5j0$ASscR!Un@z~I$Y@X?~F@zE974at%bpN36pA-^X49cH!;&CJ+~0~mrWW5y()wJHcv>=tPy(+Er7edCZQ!X(r z#-Im%q2}>U`SQnY7(DJPl;xYt54x|z{+YFKWAJ8Djx)vJ$`)`Xt`Au(UX4NT=fNdc z8|q{=4!i&91LyCTkZYe5gm;Ee@HK_nEzv{2lX_6n?i!g~szE7>psL?ZYPO&Z0~+5$ zU2Km0sB|{=D)fPOE^Xxb<1DfNJuRs7o{AP3cEaxbRCI4+kLK%o3A^AA=o~y8TNYmv zc6Z9Ly|FH~Y*r}jYn;*UT^yRcekj_~MzpeeiiQ(Xg+tUUY~|SsbVOJ&}N(}d|f$7v~8P* zZQoslFKed@hlpY5q_GzoZUqWQlQrn_c?K%Bs1r_G*P@p#zn`RD0B2Sr99d$B{ukpQ z%{>VA$7Z0b-ALFqD*)C%P(Z6cm9Wdc5ms#ZgL*Mdu>Smh*l_MWG*+*NRX5H zv2=rVbLT;B1|Q$^YGI4pI*Ju~WuSSxUIH;d47N&!>Wlb4N(w z4U$-)(1Gras}^mn?oq4!@6f)|O_=KKCx^s*u*>`;^!9cpwm#%6>!k%rIfG&z?H>4~9dlrwePP zqtxNJA?z%EA=)%;rm*5VNE-iJG_U-Q1GVd5rtKVRnsg7l_q_z+SI<(Tx;lE#Yy#`* zhg84Z7hRPyr2C<%^x&R7+K(%fQY^bt@i`SVE_x-EFE^%JzSDU>djhP7*3!+YDbS=c z3`W-*(S6=`)Z0FXb?-gs!LOU}FsulUA9A6xE76d@z#oeBeo}Rk5}f;H1|M%l)5it% zaP;hLRClox%KIWA%dHohD5{E93H_iltO_l5h0?ZwA;PHq0=8LQM%&d-3-fN7XzA68 zw$B|UEI*7wt0-04wsD}a?;)Yh;Tqb!_`Yzg?uhmkMC+4Ah>rQY(ectUT8kN?v&%^A z__dr?KHeTQ9o4-GmMQwMwtn} zh_dLJ2x(L+5UHL4vAOZC7&uQjkQU;vHJ9)k{%JS#h1rKq3+)ELnl z95+`}_|M7k%J8$)Kedo1TSY+8Ko@ye%uAZ)Gz{)JJ}1cRPjl-O;PaVhw5FFQtuVcg z`gfgZgT;H=tTqc>Zsk&R$NRLW@Fe!T(TS{v&!zoC^l@V8OZn=!aEgo!#|@)*kau#K z*f^++Cj2UwTO=VpSxqpJYm9nK~lDNr5(*#?s3h^TmQQLA3GR82ak1DB^}U z(&ht^^k<5_m{b`^yTMQ>x4t2wwDahQ;XI+XbetHTyqivSej{|WSBL>)e$wfbk;24j zypWSZ{x-wYU?p};fVP8_hKGeo=%r02E{7+7Qu{Y4;E zWEvu<{04WQom6J;fSyNsfl0bVcfK#c&a1PfZ!sy9H~u0zeeEKZds@=5TnSxqg7oxH z8SRHm^!k_rs#jM~`uq@-ZWn;#hYz$neFS!2aS8g5(4xJ)cVdug7EGU2NT;v7!H^|G zVTS2cD&O7*2R%uEVP~`HR6_<@EWHhzZb@{!QW@LHZg4P2jV}1VLATltaA9B$l^%-6 zuHO*O?^;3i$4s$P`*m=nb_l&&nt?X@OJK7~8CB=^!IqJGVD5AWdR)F5HH((R(#kZt zqOcm(w_bphBd6)?{Y*5v+Xm8GNOaP6H+=4VUN>sg3;xs@RIvm?{2QdIwx4YzA` z$y%im{!|u2X_qu=Q+*cRXm^6r>)zBN#|0j_dqMGm19Ic(c(|+C67tKFxp{scGgb$;h8;B%7dZs^T%_!=*{Q6O)ZCGvwLZB$`FaXZf2wLw|jJUSB><+ zCLHbi7t^V27o`l15OnyxfKJ@oD79-BkL_Q*q3m8meE)vIxzZ+Q4B? zzDD#k8rv9LhaCH_aQSOIuSnVYIwt#%1g^ z>7Z~3eJqdE^2cr}kA?G6XW7tU5V~7F6I~A^`9@Ss;yLWM@UeR;`NnKS+kJ0@@91<% z9+Qr(^`;2#!Lz~e_Bb@N{UN#(cY|K9&!Ca;6fRy9VezzZG^jZ(+NW)T9fo$OKSfvA z?&kf-%?GGCq)1r6JSbA>i;ANHMYD+m;Q7ibVT&JMCn?)wnP?cE1e&O8xX>h^GPTnsAe<_onO4?&)K5E`O4 ziDo8mK->?5kHtnpv3msMy2QgTbv@C{SP{--?!;!xdkEEjkKv4N5UMMl6WYo<;hd#C zYLzSz29aOk!uzSHt9(?LeV7e-xtmb8^%`N7V+i+${zgrE1<|8*6KJP}!?&IvM9-4F z(p`<^@NspJ@E@EgZSe7i7yBKAk6VM}7t#kRD(?u_yXbeLl|MY(@kKbACHUEYYzC#% z^MvJ~zTI|qB6#8(FHGB%`DyhL@Ty;?(5qMQJM!8Tz73ipw2C~XfR8Dt)H7OWPa7?5 zdwUAi8%GF(nTsU~c!B!LI>MsEU8#2A4K&=$^Lyh6P_jIQhOh4mPwy0Ho~?j7b@|w3 zh6g&HzfOi7YEbi(3AVJd%}dW*4p%)s!1-J(P2elk0oD0Eil38`hdT#}V6w3hjc z(9ZAWVcA(iwdRzVlu#?Td-q-_k9Z^&x6RAbI&_SF_4_7PXl{`Dj0vV!L*|L~^?jts zv<$j+B}?p}ddbx66dkUJ5^1i3E(dn#LJOLwh@-aa{c@UFYV_&VulqntzHi3E*;zwD zK|z=O{ZM2|tQk{g&6x^QWok^F{r!GTrp2_G4*T2xOONTZ7R-PdG9zZpOqeM%V=b9E zvtX^*-)?KxhFLNzX3cDvEwf|xtSxh3?U*BL&pNQbeJ9qDbz;uUg>`1G%#FFTF03o_ zV4lp2d9%NLALh%tF+V0TV95Mg0PD_ru%0ZC^3R~H@H%V!1b zD!axC*>zULZm^r|7Q4;vuwr(X{q5gl_gM)mWe?axR>sQNBleg*VHNBtt7KK|Z~qy4 z&Z=1rt7R|POIF8TvDd7gyxuNi2pZ)! zo64rK>1+mzV>8*TfBo5<=dihK9-Ge=u!U?9Tg>9w5|+T0vSlpsUw=90B({RBWUE** zTg}$6wQLu=(`nQdWP**3PF?O;3EF1DNPVQDO#?PdG^^)oo{X9w6p zc8DEjN7zx8$&RrscARCi6YS)_{wdC<*%@}0ontxdJiEYh*+q7VU1oVqX7sNwIA3A; ztbko**H|ID&WhL#c9Y#=x7i(5{I7qP^F4N-m9SFwfIVboteib!kJ%Ge!Je|pfBh=X z&)9QT&1zUJd%<3^I`)daX7%h1d&}Pa>%Zsxfqi73*k|^IePs>o8~e_Fu%GM~`^_5v z^_w{V`QP}Vq{xp-tQk{g&6x^QWok^FX)sNu^}oK7HfJ5C%k-E&YrzbdAv0pe%!HXT zGuD!s|La?DZpB)&Hq4S)F>7YSY?&RiXKk4SYsVb__1kmqz?@h|)`>YY7uK1%GB@VV zy0EUygL(ezdvW$=KFpVOV}49xz>xW~0M?!LU_DvjzkV;yy;%_J!}_v*ESQC`{%imn z$U@m57RCnu>kr{PlnrCU*$6h0jbfu&IE!GBYz&KHV_EdS{y5I#*#tI`O=2-@GMmC; z*;F=-O=mM$9Gm&CKa2BhHiyk+^Vod0fGuQ;*kTsXmaqi2lr8($PvpFuC9xH3C0oUk z*=n|it!3+23QJ|{*@l1pjhr{J&1?(X%C@oXYzN!PcCp=T4@+a|Z12DRKF%3zKRdt< zvP0}JJHn2#Om>WAvEwY8o%q*3$@vsJ&Cam1>>SHs=h+38%Pz7@>@v$^^1nWD7VHYk zX9eskyT%IHbymb~u$$}_yUp(W>lbsr%kHuJtb~=a2kap$W994-d(57&3ikA0zmjtm zd&Zu#YF5K)*$eiP)v;IXHLGWD*xUcz$3Ob-E&cz`z5M>ceS!y!9_SB!V8B4j8Dk=1 pBBsx +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "vtkPolyData.h" +#include +#include "vtkPolyDataWriter.h" +#include "vtkPointData.h" +#include "vtkCellData.h" +#include "vtkIdList.h" +#include "vtkWindowedSincPolyDataFilter.h" +#include "vtkPoints.h" +#include "vtkPolyDataReader.h" +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "itksys/SystemTools.hxx" + +#ifdef _OPENMP +#include +#endif + +using namespace std; +using namespace Eigen; + +int main(int argc, char** argv) +{ + if (argc < 6 ) + { + std::cerr << "Usage: " << argv[0] << " FiberBundle dimension lambdaW Lambda_Start Lambda_End" << std::endl; + return -1; + } + + // Code to see how many threads there in the PC + int nThreads, tid; + #pragma omp parallel private(tid) + { + tid = omp_get_thread_num(); + if (tid == 0) + { + nThreads = omp_get_num_threads(); + printf("Total number of threads: %d\n", nThreads); + } + } + + omp_set_num_threads(nThreads); // Code uses all threads available + + // Reading Parameters + char* Bundle_name = argv[1]; + int Dimension = atoi(argv[2]); + double lambdaW = atof(argv[3]); + double lambdaA = atof(argv[4]); + double lambdaB = atof(argv[5]); + + cout << "Bundle to analyse: " << Bundle_name << endl; + cout << "Lambda geometry: " << lambdaW << "\nLambda Start (basal): " << lambdaA << "\nLambda End (cortical): " << lambdaB << endl; + + //Creating polydata + vtkSmartPointer reader = vtkSmartPointer::New(); + reader->SetFileName(Bundle_name); + reader->Update(); + vtkSmartPointer polyData = reader->GetOutput(); + + // Initialisation variables + vector > > links; + int NumFibers; + int NumberOfPoints; + vtkIdType *indices; + vtkIdType numberOfPoints; + unsigned int lineCount; + vtkCellArray* Lines; + VectorXi NumberPointsPerFiber; + VectorXf WeightFibers; + MatrixXf Points; + vector< MatrixXf > ListPointsFibers; + vector< MatrixXf > Centers; + vector< MatrixXf > Tangents; + MatrixXf First; + MatrixXf Last; + MatrixXf c1; // Matrix of double + MatrixXf t1; + VectorXf f1; // Vector of double + VectorXf l1; + MatrixXf c2; + MatrixXf t2; + VectorXf f2; + VectorXf l2; + RowVectorXf point; + RowVectorXf p0; + RowVectorXf p1; + + + // Reading the polydata + NumFibers = polyData->GetNumberOfLines(); + cout << "Number fibers: " << NumFibers << endl; + + links.resize(NumFibers); + + NumberOfPoints= polyData->GetNumberOfPoints(); + cout << "Number points: " << NumberOfPoints << endl; + + Lines = polyData->GetLines(); + + NumberPointsPerFiber.resize(NumFibers); + WeightFibers.resize(NumFibers); + + lineCount = 0; + for (Lines->InitTraversal(); Lines->GetNextCell(numberOfPoints, indices); lineCount++) + { + NumberPointsPerFiber(lineCount)=numberOfPoints; + } + //cout << "NumberPointsPerFiber: " << NumberPointsPerFiber << endl; + + if( NumberPointsPerFiber.sum() != NumberOfPoints ) + throw runtime_error("Total Number of Points Mismatched!"); + + Points.resize(NumberOfPoints,Dimension); + + for (unsigned int i = 0; i < NumberOfPoints; i++) + { + double p[3]; + polyData->GetPoint(i, p); + for (int dim = 0; dim < Dimension; dim++) + { + float pf = (float)(p[dim]); + Points(i, dim) = pf; + } + } + + // List Points per Fiber + unsigned int temp = 0; + ListPointsFibers.resize(NumFibers); + + for (unsigned int lineCount = 0; lineCount1e-7) + { + //res_f = (f1-f2).transpose()*(f1-f2); + //res_l = (l1-l2).transpose()*(l1-l2); + //norm2_f = norm2 * exp(-res_f/(lambdaA*lambdaA) - res_l/(lambdaB*lambdaB)); + norm2_f = norm2 ; + + // If the norm is smaller than 1e-7, it writes 0 + if (abs(norm2_f)>1e-7) + { + + if (i==j) + { + Diagonal[i]=norm2_f; + } + else + { + #pragma omp critical // This is important, otherwise there is an error of Segmentation Fault + { + links[i].push_back(make_pair(j,norm2_f)); + links[j].push_back(make_pair(i,norm2_f)); + } + } + } + } + + } // end for j + + if (remainder(i,10000)==0) + cout << "Iter " << i << endl; + + } // end for i + +// Diagonal Element + + ofstream Diag; + Diag.open("graph.diag", fstream::out | fstream::binary); + float element = 0;; + for(unsigned int i=0; i reader = vtkSmartPointer::New(); @@ -138,8 +140,6 @@ int main(int argc, char** argv) vector< MatrixXf > ListPointsFibers; vector< MatrixXf > Centers; vector< MatrixXf > Tangents; - vector< MatrixXf > ListMSMperPointFibers; - vector< MatrixXf > MSMeasure; MatrixXf First; MatrixXf Last; MatrixXf c1; // Matrix of double @@ -155,9 +155,10 @@ int main(int argc, char** argv) RowVectorXf point; RowVectorXf p0; RowVectorXf p1; + VectorXf MSMpoint; - vector MSMdataArray; - + vector< MatrixXf > ListMSMperPointFibers; + vector< MatrixXf > MSMeasure; // Reading the polydata NumFibers = polyData->GetNumberOfLines(); @@ -200,19 +201,22 @@ int main(int argc, char** argv) + float MSMdataArray[NumberOfPoints]; //Reading the Micro-strcutre measure file ifstream myfile; - myfile.open(MSMfile_name, ios::in | ios::binary | ios::trunc); + myfile.open(MSMfile_name, fstream::in | fstream::binary); - MSMdataArray.resize(NumberOfPoints); - myfile.read(MSMdataArray.data(), NumberOfPoints * sizeof(float)); + myfile.read((char*)(&MSMdataArray), sizeof(MSMdataArray)); myfile.close(); + cout << "\n List points per fiber " << endl ; + // List Points per Fiber unsigned int temp = 0; ListPointsFibers.resize(NumFibers); ListMSMperPointFibers.resize(NumFibers); - + float temp_MSM_point; + MSMpoint.resize(1); for (unsigned int lineCount = 0; lineCount 0) + { + float res_msm = (m2(q) - m1(p)); + // divide by the tangent weights and also use difference between Micro-Structure-Measure + norm2 = norm2+res_tang*res_tang*exp(-res_center/(lambdaW*lambdaW) -(res_msm*res_msm)/(lambdaB*lambdaB) )/(tLen1*tLen2); + } + else + { + norm2 = norm2+res_tang*res_tang*exp(-res_center/(lambdaW*lambdaW))/(tLen1*tLen2); // dividing by the tangent lengths + } } } // Computation of the other two kernels, only if the norm of usual currents // is greater than 1e-7, otherwise it writes 0 float norm2_f = 0; - float res_f; - float res_l; if (abs(norm2)>1e-7) { - float res_msm = (m2 - m1)*(m2 - m1); - float res_f1f2 = (f1-f2).transpose()*(f1-f2); - float res_l1l2 = (l1-l2).transpose()*(l1-l2); - float res_f1l2 = (f1-l2).transpose()*(f1-l2); - float res_f2l1 = (l1-f2).transpose()*(l1-f2); - - + if(flag_EndPoint > 0 ) + { + float res_f; + float res_l; + float res_f1f2 = (f1-f2).transpose()*(f1-f2); + float res_l1l2 = (l1-l2).transpose()*(l1-l2); + float res_f1l2 = (f1-l2).transpose()*(f1-l2); + float res_f2l1 = (l1-f2).transpose()*(l1-f2); - // res_f = min(res_f1f2, res_l1l2, res_f1l2 res_f2l1) - res_f = (res_f1f2 + res_l1l2) / 2.0 ; + // res_f = min(res_f1f2, res_l1l2, res_f1l2 res_f2l1) + res_f = (res_f1f2 + res_l1l2) / 2.0 ; - norm2_f = norm2 * exp(-res_f/(lambdaA*lambdaA) - res_msm/(lambdaB*lambdaB)); - norm2_f = norm2 ; + norm2_f = norm2 * exp(-res_f/(lambdaA*lambdaA)); + } + else + { + norm2_f = norm2 ; + } // If the norm is smaller than 1e-7, it writes 0 if (abs(norm2_f)>1e-7)