From 530de400fea64634ad7ea186e9893dfe18264df6 Mon Sep 17 00:00:00 2001 From: "jtschw@umich.edu" Date: Mon, 2 Dec 2019 17:59:51 -0500 Subject: [PATCH] Start adding MPI library. --- .gitignore | 3 +- 3D/ART.py | 2 +- 3D/Utils/Makefile | 13 +- 3D/Utils/ctvlib.cpp | 18 +- 3D/Utils/ctvlib.cpython-37m-darwin.so | Bin 293528 -> 297588 bytes 3D/Utils/mpi_ctvlib.cpp | 574 ++++++++++++++++++++++++++ 3D/Utils/mpi_ctvlib.hpp | 85 ++++ 3D/exp_ASD_tv.py | 8 +- 3D/sim_tv_mpi.py | 157 +++++++ 9 files changed, 841 insertions(+), 19 deletions(-) create mode 100644 3D/Utils/mpi_ctvlib.cpp create mode 100644 3D/Utils/mpi_ctvlib.hpp create mode 100644 3D/sim_tv_mpi.py diff --git a/.gitignore b/.gitignore index 03977e9..50ef914 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,5 @@ **/__pycache__ **/Results *.pyc -3D/Tilt_Series/* \ No newline at end of file +3D/Tilt_Series/* +*.so \ No newline at end of file diff --git a/3D/ART.py b/3D/ART.py index 6dd65a4..432f4a1 100644 --- a/3D/ART.py +++ b/3D/ART.py @@ -41,7 +41,7 @@ tiltSeries = None # Generate Tilt Angles. -tiltAngles = np.load('Tilt_Series/'+file_name+'_tiltAngles.npy') +tiltAngles = np.load('../Tilt_Series/'+file_name+'_tiltAngles.npy') # Generate measurement matrix A = parallelRay(Nray, tiltAngles) diff --git a/3D/Utils/Makefile b/3D/Utils/Makefile index e56762f..e48317a 100644 --- a/3D/Utils/Makefile +++ b/3D/Utils/Makefile @@ -1,10 +1,13 @@ -CXX = g++-8 +CXX = g++-8 +MPXX = mpic++ CXXFLAGS = -O3 -Wno-div-by-zero -fopenmp -shared -std=c++11 -undefined dynamic_lookup -ffast-math -march=native EIGEN = -I /opt/local/include/eigen3 PYBIND11 = `python3 -m pybind11 --includes` -ctvlib_config = ctvlib`python3-config --extension-suffix` +PYCONFIG = ctvlib`python3-config --extension-suffix` +MPICONFIG = mpi_ctvlib`python3-config --extension-suffix` -all: ctvlib - ctvlib: ctvlib.cpp ctvlib.hpp - $(CXX) $(CXXFLAGS) $(EIGEN) $(PYBIND11) ctvlib.cpp -o $(ctvlib_config) + $(CXX) $(CXXFLAGS) $(EIGEN) $(PYBIND) $(PYBIND11) ctvlib.cpp -o $(PYCONFIG) + +mpiCtvlib: mpi_ctvlib.cpp mpi_ctvlib.hpp + $(MPXX) $(CXXFLAGS) $(EIGEN) $(PYBIND) $(PYBIND11) mpi_ctvlib.cpp -o $(MPICONFIG) diff --git a/3D/Utils/ctvlib.cpp b/3D/Utils/ctvlib.cpp index c5d0cc1..3406516 100644 --- a/3D/Utils/ctvlib.cpp +++ b/3D/Utils/ctvlib.cpp @@ -46,9 +46,9 @@ ctvlib::ctvlib(int Ns, int Nray, int Nproj) // Initialize the 3D Matrices as zeros. for (int i=0; i < Nslice; i++) { - recon[i] = Mat::Zero(Ny, Ny); - temp_recon[i] = Mat::Zero(Ny, Ny); - tv_recon[i] = Mat::Zero(Ny,Ny); + recon[i] = Mat::Zero(Ny, Nz); + temp_recon[i] = Mat::Zero(Ny, Nz); + tv_recon[i] = Mat::Zero(Ny,Nz); } } @@ -77,7 +77,7 @@ void ctvlib::create_projections() { b(s,i) = A.row(i).dot(vec_recon); } - mat_slice.resize(Ny,Ny); + mat_slice.resize(Ny,Nz); } } @@ -121,7 +121,7 @@ void ctvlib::ART(double beta, int dyn_ind) vec_recon += A.row(j).transpose() * a * beta; } mat_slice = vec_recon; - mat_slice.resize(Ny, Ny); + mat_slice.resize(Ny, Nz); } } @@ -150,7 +150,7 @@ void ctvlib::sART(double beta, int dyn_ind) vec_recon += A.row(j).transpose() * a * beta; } mat_slice = vec_recon; - mat_slice.resize(Ny, Ny); + mat_slice.resize(Ny, Nz); } } @@ -183,7 +183,7 @@ void ctvlib::SIRT(double beta, int dyn_ind) VectorXf vec_recon = mat_slice; vec_recon += A.transpose() * ( b.row(s).transpose() - A * vec_recon ) * beta; mat_slice = vec_recon; - mat_slice.resize(Ny, Ny); + mat_slice.resize(Ny, Nz); } } @@ -256,7 +256,7 @@ void ctvlib::forwardProjection(int dyn_ind) { g(s,i) = A.row(i).dot(vec_recon); } - mat_slice.resize(Ny,Ny); + mat_slice.resize(Ny,Nz); } } @@ -427,7 +427,7 @@ void ctvlib::restart_recon() { for (int s = 0; s < Nslice; s++) { - recon[s] = Mat::Zero(Ny,Ny); + recon[s] = Mat::Zero(Ny,Nz); } } diff --git a/3D/Utils/ctvlib.cpython-37m-darwin.so b/3D/Utils/ctvlib.cpython-37m-darwin.so index 8e074c34a310ae5ce580269522ef3d8a41bfb9ad..001a5f41c6b1bdaa35e931aa86087d69da659b0e 100755 GIT binary patch delta 20760 zcmd^{d3;S*7x&NJCs$%hLS`|<9Akr0st7S8Vu)EN8E*zM-Xz3O_nKqAT==a2V~_k2EQo!@t@wTH9MKEu8DB>DfX zcYketuQQK-d*-ZZOxHB62CVg9_1)ujk@xM+G%YHsf8Wu4qe6xc8#NjNnpVc;Msru2 ze>2uny}CQz(EU9B@3Mi+F+mS>G;{K7_&gYRLseW|Z{Q8bNT-*ZJWmL&PQdD6L3VQV zb@6u;nw)y&{oPj;m9~=l8ce=_adkhP9AD|JTWS;yz(K$Tn3|FZW299$_Bpll3V5LE z%?^?Q4Oy2xaW@5ctZk)AL5Ji546#OzXj6b!0mZPlxFDNoAN+iSsJUp05@$PZn8uF> zh=39iuq|`Zw4ffk_Q+Gy0*C3^0XH}-lXNWu#twVt3+wm{UF!t(HGE&!#+PYYZ7d5b zE-fGl*AUQ&j@31qn3A^X+AeX|U-xL`4^FJ0jeAUhQ#l5cnKD0CV0$`PW zb1S-2<5lJnk)~@dP*EH3L)+MBQD2u;&W;EdSEtZrrnepE-8wQ$k?#JP&SpCfdo*D~ z95*~7n7^a5XH37O);0s;y{ltpF-duBQ60AWy zhVuJT{)n#~-+FFkUpR8Tdaz}Vyg*Y3ySu6 zdotHeB~W*d2$$?e(8ujX7k%6rGZ&TnPGl^ssDFbk?Y%uVLdLSLe{F!x{3Gt6OWZ@3 zd7UnLIB*wtaThYy$?-v0s7vvB2{(%V3j0oHK1COXhB>qU76pvy&iG%;923TF<+i0o z&&Fp<|4|LMMSgIMLzjiXZkxUaZ|Ji5cqITjE-AzpUf2CA_8@m6cP&ZB^Agp$v}*vLPZma;G}@eF^YnnBn?+o zs$yZ1X1DMSU09WE3axU>ZW~%PRLiLqx@>^M?Pe{HZT-1c-ad@ML04n5vA?7J&F*@1 z6>}us4B{>O7j3v1>!g46m=!(y*@N+-KZ`ud1>+4K6@5@qi*eV7jsv&mv3`!8xBJ61 zu-&fBc#r=&R{XjYCYbMUqZp68=SaO%$`4gK=HB)$8uNP(#=G4uvRB6Wu?^yWP1c96 zzG!r<#h~DV3r4@%3?@>o3u1O1)`|ae-uSQ%JIi?gZ^S%Tmf`>YS*E!yIEJI%v)4ii z#}(K_T{rd%x1SRa+}IY*&Km37Ss3Tr&xo@6Y%}*gZOrjtO*kJ{U%chTCi2JijK92C z1mip0#V{XM$}hNy0AE(j+KU^$Y!zSNDlC3%0Pp20KK5gs`FF>~Z+@&VUw>Q#HDGP{ z)Z@m)28?kQU`%SroEeY!!boVu?lF&cM`1jd=%FjT+WhxfRp^SYZI2i!O;`bA(MH{- zEQIM3P23pOjAb%CU+2d0=Ik}b4<9s4Em^n|&pu#`Y{PcxY{mKxtOJWOI(A^0CjQ$l z0o0xy5~n+}0d?NK&v3%$z>(A(O<(d5ox3osX*~Hw?g;(4bbV0{{(~4V0vtG=Hu`9-_ zz1f?LwKu-&!vb}_^G8DuVeyO){6Qr5XW{(Mn__=|R>cyH{R3D;UEXcA;W7#iBpWGu zj$t!+?n-0-81@choyEj)z@f{H3*(UMi&>GtzDtb~V{_((pR%3lUGPIun{AjkZKLHuqPXaEP zX|Od$!Rn9U0os0Cuq$V>6MVeGBlVbYof; zGPIuqyd%|EpMwnTPY2GOYV4nm4DHVV4ox;L%s__rX94>r8V_b6L;G(4m&Y64bC98Z zEAW9>qniyG+Mf%Y7j2B4$MPAwZu~xKAR8Y3=$!CrIS9fK!bN{a{xZj3Qne97KoQ}0H>sW&@LoqF!i{yz2U!LttJ zor5r*dbeo%2W;WgJBpj*refI8srMO3b?RYEo_g(}Ryg%8z9LUOdz9!{!dgy^fpGb+ zv8us@zy5A;F8HCrez3s?C(t$)-AhbeaC6)g1RH9w14z{%#$(TpRg;8JroH~VHP=zPoA=KjB7VU<_}PGYZaTx zGs+~a>&$2JvhS5X!k4hVVhH1_*dB2eH%*Nc&MP>3#TaYiE+*DdMAhQ&un~q+ZC%MUBS%3~Vkm=5fr!c(n=N!q`6JQUG^j z{N2xGv&Pe=d_H$+e-H+DS=XT=uLTd}o}U?eTk!3Sof46)cqp4MwzlG98g8i4wDJO& zf}zWLef+UJ8h>1RiuV?4mGEfIn_Isu(X^_TCDNXd_LQ`zr7e~AoV4eqy&&x+X)jBA zRod&)mPvb4+8?F;S=w@GZ%O;Bw0ESfw93sN(%zT$fwYgLeJt%?(ms<``wF|urFD|F zhO{oy){?f4w64;+OIu%BPieiS^_8~4S5~a9p@hcL21wgX+7{Bbk~UD-qHq3+fUm5U%}H{Rm*`ALZuxn?GR~)NgE;UNNGn)J6789 z(oT@}b!jI{8!c_DwDHm=NSh?>6lqhWO#^GKYMCLyB5k&`)1=LncBZtmrJWQz`8q(%zT$fwYgLeJt%?(ms<`J0Z^k zX`Q64A+3wFwWO^ht*f-|($<&OQ(A9neNV{w-#|h`X&Xx$AZ;^gTS(hV+CXXBO50xA zj?xB6+gaL|rR^qdcWGah_BCmHOB*b0zY}u)_m?nG+E8f+OFIN?QX1S!Dn`N;ENKy{ zP*iWC>W%6%R4=1Ci>fWETc`q1)qwj@g*U1usOq5ViHf6&KxKWzL8hR*gUX8PM^u|o zT}E{PRVk{IsJ=w?J*tDK@DaCSFDhp+74Qj7OKOGcO;oRthP~oS%q{SS@ zzhWB7E$B=^wGUM+suQTjqxuEaFjUy}3i$M;B{f2Y^(Vc8DhSnhRM>h_CaQ+0mZPeV z3clgcDqK(c6O(p{fmEQC0Y&YJtiXRR}64RIk^7`B(86=1Wp0IxA5v zMfEeP?WnGzI)UmOsvl9ook&Z1hU#-ve(;4?#Yd>RqIw_IFjNJo5>c&3Wkt0dl>^m0 z>D!4i7iBRjobE~Ap^8U!2h{{r?k=G4d#R+CPz^*i1l4P(l2PFoPf3eV1)_QjRby29 zP{CJVa3E0CLiGR@9DjH^UC}5LK9uzxRn@X$kEYp^TGWKl&+w@)DFM}0RC%b*q52S2 z398yqQN`z|UPkp1s&S~^N3{!80jg3|>ru6WFV`xTqZ)*2UR9N~D#?m67oGc1rJ?!- zRXi%^I-n+?YK1BS)gV*@Q6-^z4b=yzI-@#^DiGB}RE<$Js|(5#)o@g`P%T6S_g|22 zp?m-j+e!OT{fg=|s+*{8qq>C3Wp>Mc|q zQSC(40@VRj4Nx6Ng`X*szCl$3RA|yolz%-1=Ur5Ppkml?IVyKl*HPu5I*)1rsuQRj zsESc-HKugpf3|14jQGj?Yn}h@5bf-|u^17{AMuB4jPN)GvBmKv^&Q8 ziLNPpGAW{@lP8wp);Uj{OotXOFBakZd22B^gB!f_Vs!Qr!I|8T zn-+_Snb4QRi^R4}K8>$jB>XK8=xK$XC&SjpQ$>Ic>iTJ>0kEC(zPZK^^SCSH{ichli(nMHFMv_pJxwfI z02McwCYmpVz5Oyr%v#9r@XOf-B#q^KkwqNM1Ca%v{+IA>e27JiS_0g|B6cp}*ZGu8 zW7$&vfblKqV*Ya2+n989JWMO#cuYzc6IVd>k?A-tjYYkcygP4_F2YxWa~e6}OS0~K zY`W;Vim&8%)5OVDu!kR};e;8>=cS>)mngQw&WuSDSL`s{EoHEmu&oB?qg1hJHSG4e zRN=7(YS;~&$o<5@HGC4!OBMbO{w5!kDn50Dm?ij^yvbAy&K1YlCcHvv_v{PvH-f#mQ8F z0UP0X7|F)6jU2vGNK6vVH^V-3ND{L)!`RhG5^!9biElSU*3(4czXiH|B2nx-!Fz~P zTli4E8=P-~b5vDtdXV$}hv57=PSh?0=i)e(*aJa*XfiLnB4e7cc>bM~d5Dzzl+I_?K|`N{O2BZj(%cYDV*=`MHsi? z8*s*bh8suE@!pK{VPfKW-dFTL4~0J+iZGUE%WW^C@;pyw{93q}cmW#R6)ukAR)mXo z7h#vChGYBf#GQ+L8J`s{W?cejAUMB+nX&desPEgs2yH_FjS*QrV22K%b9KD8}o;L{3jCMxH z@A*uoy9DCzkS+DH8yvPBYP$!*BQh@1-^im9AFyEwXsxy&#l&BDOCEef%=?9R?3`5n z1FTGgpA+#J9+rRmZGZd4&0l!$<`MrYr$zqVKy`WS5VzeBmU1|p`8UK?*jh`dTy!IN zs)ePnvD!)i<8t$#^}vo*hvfEuWUKZ+R`5TThi9n&rycNyLHu_IDsOQ418-4;OGhUo z=mrm{VPxFm!K`-mB!Gh^=KRhZiZi#lyYb6y-lnBq2s3Uj%=3~PM%#M&ehxz1^w5W) zqqC=;hQAlP8wLJ)v=dYu>{O>-(;X}$AR};O$jE-eJ3wz9Frt|jJvAiSZ*X+I{?Uzt ziuuT((!m+-8HJ{d(2Pf}qr~t49X=2lNdbBh<91y%YNm(kmP<86Ml=1GQ)VHvi>6_^ zw;0=8AE8?gxES_sdIcAbR(gG&Su8$mrT5V-rIF%pEB#gH%uCE}DRmV+TI*3YGs{?k zoqcK#7Kc958#K0*xY~oGq9%+E8|Gnh(O&_ zcz>>YaT}bMf%;fq8$6iGg+*(!iujm_X`00@+IaJzyL8`{cFp2aKR9a8@Szb=De37MQSmdGGSWK#{V7yq- zM{nR`iT+2zF!OZmL6RtjpqC{>2Ky(+Wk-diXUxD*uDA;!-C|?HHCS&j(vruXCsA?x zKU}IK`7ijYccMV~23RanA+YxrtbYd+GlTVpR!iaY;6QVBToMKkyx>c>q{hfyEPn1B zhRzfiq5-*a=8Wv*^fWA0@;qU9TwIRDVva|5>2vo;b5?e`MUKlQcfDY8TPF9rr$(o>&@3b?51tNWWQJ0M4N z+M|1zpxHt$Ug@VdhqIH5N&R%+=9XerMObP^x+ObGZHA`D!)%emgN0Y_(Y>uuxKyPM zNzQ^^WMr7rCl(I;Ho~%mE=oLON2ce0vkCx$~)AY8Wp0> zeQe*R%V9w8K#N)SFOOWZ!%-nHhS}!MSc6>!hUKKj!U+_TWR9DP?)-na<$$8Q;2-Yl zUhGib;jvT9ah@>Cp(j$ynNWengGV$xt)G?j)Y}*Zp)$zLW49pf}Ym z1x}*Z0NtmNWrr%Grt+w4I22KI;qr$RA<-}%#!hL@9LZJ3Y zI349h&1MoY1N9~?EO{n5PGi!N~JFk8or9@zc)A+$qS9QsK2#tD&c5{CxjrLe#x zu7Vh0*`Z30GsjG|m=kbdMw%0_)j}1(T^wt+WYKl$0D9pLKt~0Fc%_X{2d5aD>b_5J z)DX&-Q08ddpCRe#Q*$z~+ESCS4bnYzOJxlq2I+lFmOSTr!J=$GTobT^&SFle?jzb5 zdTrza;vS)TKi#rJQY+n3=qv)lbYF2ORCj@+M|925dtjkr$pN~hM2ZGtQkdQYZt;?Q zL}8fj49^G7;*&7Ft8OWiB2d)%Sg(nt%ZVcfW179Dm@rsxWU}Pfgl={^s7|bcnqtFX zy)pEtruamP!kXgBV7;$yDXuA+gzL>rmJ%#*>7d>OD=n=lriAN(R!cc};SmBJAhOMt z)Z{dnjnm|C*YW^=R7_Tu8LoGhYj#-{RAJ}UDy(boPc9JDs-m~+bd;j*?Ym3W6^j20&er*_$ICYuD znSf}{m6K=(>IAyh79yl1#-ygk;6NUbm@co^#i-%#ice_}Dv@YTjZ=>WrG%;J)6@lx zI1A7C&Ty+hPJ@S%m;^XhSV4JhcqU3u&A^KS)SDXxy{^7JR@N4KhU$JM%Ogk<0mZs! zXDDH-19fC&TjEk>$MXo|)yZHdOfk#DonOcB8m9l@R-^1Fd|S-LnO=HV(QmTe$JJI^ zM{{z@lh$5pteUJh<@FY{(lq=|y?R+A4oB-v>fHi&3;4KH{bBA;aW`7;V{P7A(;73I zRtn#6;?kYCl6WX_p@)iRSLYLNuFfakOKiul;bAEyE+D=`T&&ot)eWQ$P(pj+O5$0> zg+8i62eEu56gf7HdPsNeI za*hm8{57#HO!1xSIB~yDDnDSbiccW65x+~E5ia*%)2@@FAY2uAK&%Z>91K0j4nz=d z1ja*DI#k7P5eE-bJfI8h5U~K3<%g^I3F3_5IR3Es!BAjA871@~E{IU^J=Huy@pr`T zBNaD>k;C$pqZQ90&KR#apSXsj>~2>)o7h@Oj(5R<71-CR1fL#?-8U)D2gVKtZB=}%n%`0U zJ8|@O#dTg)`Gp@SZcA+cNbv~8R;}Q$a^!+TcHorad&FhL*)SO~URtW+JK^_0WcTxm zUEtRy(gy=sN&mEw+;IZQ40~On9!bE_@XOQAzcw-Hh zcZuz_6;}~g)>SM*RemW>Bv^h5h4UXfP*zVlD&PWwZ1+^$0WKJly%dik4)#|(gE*j} z;?=~3#Gerd(f5Tv5LbVH=n`&K4OD-T*a0p;*kCDru?SZv$%Xh{1uT22;{l5ARx|z_ z56jzdfx`099TaPDfs&_nC&j&hF`m&`@$Xi0XkC<}C0wyEA)>3|Xkt6@k?J_U)B#IP zxL{#^z!=3dh=a!|zC|1zskjqd&@kUNQL*)HauiHb4*cmI6Yvd7TspxO54jB9;lyPI zu@LX z@X;z?%TXLg957w6oj7=g;)}%5vlKTOqsrUfK!)SLiW~*_iUuqfh>PbcZZuXM!?O8` z6M=CL0~RRWNbHVpd%$voSX-*N^Eg#Lg4ln&;?h+RhxWDI;J^wi(E-c+NX60kA_*+x zCnzo^?l6%$uv*3W>xyl}7m3S=KbWNA!E03h+{uayfw6sUY!p@KPzmj#72Aom7{%IJ z6~9PqBi;=+Dck{Vor*7nn-sE*_;q+3mGSi|-o;Gu^-Q*})k;tfccBt)5El?1NL2Bl zd=*~>cR8%Em^dL>al{4{-#~1~cbi~2OI-N5;a9H znKKoa5f7QAID@Hp>)DDmu9&}}*iIZS8SejOx=P5MqY^4jivJ`o#uv$8nE-QD&f>a? zV{M8H-4rLzRqT$hpuv(hPqEEYaqfJ@K|YG_A>;kOk{k~fs07VVB|KTEIJ$x2r;8NZ z{D~JUw&UAxu(V&IxUjL}cw+Ye#eJ5lctBId4VS_7#}#x?Gv$~D4&U_%!j)#P^rO@y7;hwNMrGT%j6hM%;0w;_bu<#CKb&d@HgJ zf0s!=RBj}${*3z{WM??5MyvdQRf}FT-N5o;o&SO-C>BR2D zdx-;x9}x!=*IA?LizXg!B}Xnft`p}Gk8`L3f_MdS0r3~ag~TCiRemw?ZsHQ++Ur#O zG_jSfSB^{M*h>5}@f@MzmBgle#m|Tz61$95BeG_LihC2+*r+&wcrI~}VtD@BtP=WB z!d2o>;%~R8_;}(zZz_%^&L>VHc7996Gl+u=6z36tN}Mklj=%fcDq$NXjNhubkoY`t zF>%T^6~9D$l(?KY7(Nl`?Js;+p?aycihU*J8=Bk`fL8H+F{41;p2h%ZT0HS8;dxS-%T$ z5b+4&>em6%h^t=%SVvqq9$yQF$soY+P@cNgsd zaoldjrNjXrDz1#Qs$+7I9NGlMD+^VEJMn-$iUWwNh`ST-->c%G#HmG!#}l{sNU@E$ z%Ra?rR&xACjz`4r>{kgc6IF*M8j1smeGe!OCjOQlUud4>Dh=Yj5XArM#lYe86hwL-r2ztfh4srFb{`Egs`8J9N5$6*RRczJZRRF9{ z%cg|tcRccmD=A(`T>Xy6Wn$l(s{AwJP~wJ%R0nd2dl2s+o+KIW|E1(eD^>-{h?^c( zTuFS9So=XWxbKLHyA!7!RUAOv^qArx;`7A8#EZXxxD|#pf*ifRR0X1mD~U6Rw;fk; z8}Zl@#dhMDuM`&$``wgJ8N~s_%}Nyq5u45`4kq4k7VrOcHTmxna1QK*1mxZiR0E$8XFODVhB%w}25~O& zJ>uEKrgN%38?g`Z0^&C3;QE6N=8>Z(B`hZ%N^B>NCSFUNMZA}I5pgN;2I42g9}stW zq;}|#l^m1FagI2a_$T7`h#wJuL0tPAwF4K38x!9r?nG?*Q#IJ1xH++P964Sl#}wiy z;x~xni0#CQ#M_7$67MJWd#oD#inu@VRpNBTR_%9kWKjY?ul8sYEY~sd};rK5mM>9&;NZf^Z7x62^#l$^{zaj2J{4;R~@t?$li0gc- zc5pc|w68TG#}-QHOnj4g0P#KI@x+b)QhPX+xEb*r;{C*{i3dMb`R@=X0b~1`L5^|H zRKf}3Rm9hbKOwFpzC^5FPxqy-HU1uBZ&iu zcMx|Wt|T5v%$!txk;K8oDa3hJa#+c+op=qgLHsVUw@FoafOsPDN#f(g*NN-YQ2BR> z9mGzT)DBr2II9F7a*QNyOPu4P;=PDx5JwQdK^#ZCgLoS8+*+#q65=()n-yENyX5$g z5&~e@1soJFB7TqfI`OB(HJhsZGsF#ue;{s0Y<)c@lgKx8`z@Nk=#Ld1}@#xm7!8GDR;Kh7S?sJt{0uw>OO3XxGVG~Y`dtJ zulHp8#ZcUQDQ>LQTbs^8+67T}o!*WGi5vO4M^pTh;7hQ62kT{6ufU3bXdD6SRai%e z)OEUlJ^WztGpzV$#O1J7h`t;2=4^uaY@ObE;63n6g!MkG_&|=&$oSZbPmZs{iqBX0 z%!JQLlVF_;D?Z{x!5R%~46LzSv|X=vtBV&5e^~Lt&`6}N*F);Q1ac~*WWqXK9A2+C zw>m%ySl@*8Jy;D`kHdNn)^A`v59_zE;`f=?VEqZ!|H6u&f8tT;G;f$h&1B#-kpk5F!f->C8Bo2A6V;Ys@RSuaM znVG2-iusV0;Y_K8W{N|0wY~~zMUzsg{PsR)llJ+p-*2toTHm#P|NNe{&VJsnz4sZP z>E6quzi(_-@JCahFaP@Sp@*h1UDGsoC|*zk_xfDqVP#y?%;u59rwj*I-R1%(ZXwcY z+J7*%g2j`ch_#^&ob7`XO>GfQt5Q0LLKh+C=-P6W4zTh4uYU+m5&wk||%AiFt-d-^*| z+yZ^9Kd6eCL8>CME~B2Nz~sm`vfHWz<#K%l{L|{t`l_43bVn{`9Ud+)wLipVCVkI0dANn;n%;4qZ_<-3F4nz1(b+2J z;pT1G6zBEk5v+r=w|Cr#y@#U01obqj7IluPa#?Cl%d%SQ?a{IaAZsifL zuytNxxjq+)?~Q(!d-@3(8{phCCd{)eUqV&!{V|tx*0%UU_!tlNpg3@9FyjFmoY7I+ zcwt`g!>G+2Sf=xp#XZ=~;**Pis>dc4A6V643_sP%dF9L&KBQ&wqAyByzQETx=-d|m zjd$^-b6wn6Kry?xfwR5E8!r`dUf8g>|JDAC=hQDw{iX%u=6cTOt`+d(p3ZjPzQVhC zIzRq)9e>@!`HZW3qjdwhwyHGK2NrG3K&Qjy9gql~$gq`nx>oKO%e2*{S*^lW9_tzd z2~h)`MXoseZv(hiTfVZ!y~&!EaDc;dLt)ycMo=v*cV#eo)|osch5CVXQvp1afet}s zNGnrRt*9%U77+VEkwZ~KRWDD`Ru$XHH9MvjD6Ce&bV;FM5J*!2t4X%_hN2=xuM)wn z1Tr@%3UQOPUQxMNo!`m}`~Y$1phP zYHcnDJG`Kg?rtadKH?OQzc=K+jI-7C(k zj&I4{5Wm-B!}#8BjDGbQRNQ{m7}0>iiPY_?SlW>F;vQFw_ZzaajL*F!Ry1Z-|9xkf z)@0Q&T=hXiUX*ZB!CKL%3A@FMD#YC;>^0u*qOr+~jp6*n*TU6|ZQ;S^jb+VQTh4{M zc-@D^@PY2eeIFLVc%_?|(1Mk7UpEox$4XgWaovxt=f}86Z^=gSXSw*WCF{-qJR^Q= z$%gX-XGBmd_5{y8W5l##jI*xBtN`Z0xb>8g(3<_mnveSo)^nX6_M*=de>|)Wd$Hfx zG9$SyD`YIuXw;5{GX2#*nGxBZ*_d8&j~ORBu)&O9KVq0Vv2kvE;~``66YO1`ZOZS# zda!v$&mPQX;(xp=fV#7j;)~vFRKx3)3{UthIFg#P-BZm)?>?+$y$`_~_TrFTm10UC zHf%v0VC9~B8l3)3RxPZptv&wd&!u=@QE4M!cf6*lDVp|k(7n#o2X`k;L9m0hMcD%+ zZYO0-(;UMpXF{tmLyyyPhFYcZun*ho?*GYKcz`_cL>=T4At6A5>6Gy~+&vlha*1yT zvGME|(S0y`yMf=g&~#>h#}Kb?#kIk#L%ZA8z|my&_i!JVcfs74oUqkbnyl)BjsvDM z*NpZil+vp@*`$jIX&Ul18#|e85F>Fz*M%Kv4J9H8gSwU8Bwvwo5kO8z&Xzu*>jK=3;(&mVXKY&xyaCe0mm0n)$k2Z(@btw-|1@Oi-wNC{$B42a zL;vZ(HCaY>1~T-Y1zeV4HVECx(s=uD85HZCERyXrzH`FG%sv;V^4WSEpn%;jaT|9q zyaL!gZ!~y?+{rup!%O*xzLZoINZ2ir5_1-dIw^lH8cp_^UeN#AgOv&A&cn{BE!+24~B)57|?mtclV12rJ^={ioq|)@0Qh zIQ_=2y9pOp_i4ttW6T{McTwW`Wy~qcae08CGgg*jkBlNad(z)h;yu^&i8nWeP6TP zzOUDZu5$aW#oa$2F`|u(3#^feUu|rxxWazmtVX0{Sz!~0tn$5>mD?SgB#t&x2l zZnZ2xEV{w&vQEPK9s877jG*sXWiz(aXnBuW-Pkze`M=o}#%e{v12~IbH4Z*t=NR|A zE^Oa1Z?UbGE#^5c37dNJ#k}%crH}A+Y>XJs_xg*H4M&iH|C z`D=`QVpIh3CXBy(NOo&HXvbG_&%PhS;^y|7D%N)7U3uT5#=eex2V>{OjLtla<%(^c z`P6`_-!!eV5Kh6c+#!cQlt<&&iU;^?`SLf>ybJGOzYE^lPBqg0CGCA_A4;no#3+~6 zOT+|C~dy91=4Poc8j!I!P;v(6-wA9?RIH*NV`+oUDEE6 z_C0Bfq}?a&erXR#TO#d2X+M_skhG=J9+CE#v?ru3Lu>E!xrEcwej#nSwCAKfFYN_s zE2O7N$ZmKhP2;HdsEs+$X}e0>P1^3# z_LMeA+TPOkmA1dM!O{+tcCfTVr45mGgtQ~29W8CxNAmn1D`C8}6QqrhcCxfnq>Ykx zy0p>K&Xjhxw6W68k#?@M3DPD?J73ylX;Y=ON}K+XJpVH#WJ{YP?P6({f=x_^M@iLW zxPv9;q6$N`9Tj|+*AkDR>Wk_Us&1&@mZwz(qVj_0P?ayLZm1ffdIl9oH4Bx!2EIZi zE=Ktis+UlGk7_rn%c#mwm7}_X>IABvQGJX`2UE2VRSQ(`4NXhzk7_Hbv8XnpGNXDP z6@J01LX#*stbbKD%6HM3jOsY5IjFGrs_Cfip_+iI84SH@B&r^$;7cOR2~`lPL{yzo zEk_lAY744ns1Bm?M0Fn3Lzsg-@h_CWgG}^+uYy%Spz4Flg(?Eo1yl*BPNRAj)e%&A zs7g?6NA(`6gQ&Kn`V!SCaMFduA;h(>Kv-yQI(-;3_q1s9YXa4s{N=&q1uINHmX8Y3sL2xT94`lRBxhM zfy(|7${dvCs8Uf?p_+>de@Lu~MimI(%&YK6sl;KZMx&aEYA~ursPI>(#C%j;QN53< z4XWd)yiwgiRUef*{4Q1n#~(hOzS#O~k){ogtnKv6UQNqO?2DnB@U1T~8`V`*1*p!U z`V>_eDt~CG>JX}tsP?0phiVt9Pf!)2a-qsc)gOLYt9k*|G*l~UYwfj(>rv*Q^Ej$h zRQFKLMb)AqsAyFEQAMDdhH5mb98`l*eT=F%s_Uq_qVi}2stu~%sJu~`QPoGa2^Bnl zLB5ajE_~QdJb~&Cs%xljp!yqC1uCD$puRxW3)L}H;iwLxibqw1Y9Xo}s9>qJsx7G6 zqS}CJ7^=0XX5#u+EkoHGewVDuL^T1`d{pyM#iDABiynn)EDmNoDhH|&s0y*WU{vp; z>WS(jR2@;3p=yQdA}ah&k@y`dcTi!8zoNYV0GtfF`xRAFRF$X#QC&l|8r6AJ8&Q3Z z>P=LosP-B2d+|TJvk#29v-ugF4_hmq%;Rmuq&d8XKl7Y1E}k#r{N8GD%ECwU;?<&g z0*~V6)y9GZeu44al_EEp4=|P{@>z`ETQ067@q8Y?Tr^m~M~Ra8d=mfCfpC-Gwu>i^ z@e!h5GM~by*$t?7i1SyMig9UtC||Tx_q%M2L;I ze6eIT|A}`?G9W36^9u>$=voj*6NLXd-k)zt5RvPE*CmME>-aT( z&|>61&+jt+W4u`T0?c-AygDAH7vXrk8!uvBgy!Fg$8~8Vn!dz?d0M;}_YycykQ06+ z3+6lG#h~^4B|bb}oLUbvtet}=OcXyi2mM1tX&y}GtvTXK9#7y|GB`vyUIu5_9P#SQ zFzxnpMDq>M!oxW6>O9_3e7u3r;uUei-^sUf5hp%!!m>Av6S*59sWMjhZ-V`?Sh0H( z#6OJ{-^yq3gE8XC*KrYJ#L?GbX_m)`Q#m}yh%Dq=89yH_ zetQGfy&zij*#>9bof+cq3^)b#H{s%FGDF;c6RJi;iv>Gy6}Q8*>P3qy+rfEp2DSq; zcnfAQZ-!{T1IF5R24+U_is?ALAo1G{Sb(Rei;v!h8WUs0->0FOyqz#CON_X^6P#`3 zez|bl1*>x@TKMk*r!yLZ%SGSa;0%uz<9377gPe2sfb;MSv1AW8m%tD)bK|?<>^DQ~ zeixiSPgl-k?}783>Egopo2%HzEsmzXh!TItuF=#J1$B}c&K5&|* zi3R(>`Dv8OTvZItMN#6_VsQ2$=l=cR{Bf!{wI7^sPsPkA5w1m&rd0~X2(2wCt;%|hjbAQsOJc5rP zIE?V!copeY3b2ac4TMevO%DT{7^m)8&mZQ)_^0EX#6mHVQ9} z5n@w0T<`Zr8n?^w20ky;IC_o`Wqk4o5p$jo7bDL@__5*W9AZ?T=ShtB2oW(C zpu-2l#8KSe8YZ5+2vhoK82109_~|0g<)?;;B^BVz0p}$=T`oa;oh0P*b3-xstSG$< zr?-8mxN;e;p&3K*LTDpwSNL?^aHu$S1&-p|LqyC~2u2T4!C$UIu+b1P?i+C48f+Z> z1_ym^5MCKi8a=<|i>G=h;+9;P|X9~8&u;5ZP)_4 za53mQ@5Ezoi84et?4z8`|m< z&~c-!p6X%5cF|+qe4Dl_V%9)w*U6!=3qoUCj*V?PGB$9yXc(k-@U@ol$wB2~tzOn5 zlQqm***Fr1B_;;xe$3NI4AP4k-`Pa8?yZM4a1`8wkG;4tQGb}ePqgZzH>>X`s8P;Z z(X)@9q&o`h8L^}FDlS&`*PAtRJkb| zd>|*}jP5IXKdnF0`0>!rM#BEI-nx^cz#}ijY>u8XCNgutoVdB>xa8!tc+FA%sFflg zaqwxqdwoaQD|I8Qd|H39zN=xZS78C$4}JW6M(-v@2kRaEDu1c#Jyw$q&y9=7)*L0G zB3vIUP6z8gZ5$3ajF{);VD$KU3!4wH&FkfSz|}RA9dXp7J3QKa6OgR{32S1>up8ACwf0| zI$W>M9XoLp5hhy^Zx$jZeWLq$h^&MJOS*ZoCD{^}VZlHJ7ey0vUz5!R8RFU} zddo@Rt)%otnU)N5WYS`pTCJR+qq9=uGn3L%%@ZsON9AB6H57@nW@IH>%+FXdL(@_- zEje;P4qXII)LZnot<~iiLLX&Cq$NRGIyRB_?;tdd-U3}Lo~XC#Z!1*ckxB8H=Fl|j zA`I?)%s0}KJSNkU5}BEvlsXTCMGzE+59z*kfRe{DspZmtxztF@9`n^rx?K52rKg)i zVbSDRD<1hq<-}X8xESGanemC14D8zVC_LJd*)~2Aot2L{)6!Go8a}Q6p%AJF4;Uh z-jYO$1HW{dHRn{DFP4-m0a|VUkNVJ}Coww_4%4 zzy(}u63Zv+-gaBwW8u2<>oHf|5icZPTIxKQ!+6+XD?61p+>&ArRW~4T7Cmx~PPfRJ zmXIrLE?xVan5`TUKMzDPP4=8?3>3PHiBt3d-L}(Rq)yR&``C*9 zovPAGlrz#|%d(`#Tg+3ktm?EWg@71-RQEB#n9AJ6*(rLv*0yq$u1;Zf=D>t-G8OKk z$yB|?K%49F6mw=q20b(&r}Dpr;^jfY`ql0tb}HO5Y&HKaDMR((@Q@25M~GKT-7>Cy z+)pQ0-s3)4mAc%5fBEU=T1d`CR*Sjwn2gA*Ihi1FP&*%sOooeA-usHcDI$*P{p@*wYQb2I#W}V;?w3@%%juOQbzxC*D1%k?;nGkKGd!75PPF^-f=3?Q(y;SN58;9qbz8NEco3zxgR9L$w40`Tx3W1ro5M45YU;uyc>chxDJw+| zx3Hd=I87g>+e+$*m#67%bz50Iae12VFTR+j*VAp~l4h7}74@2ii13qodmM(Vo|rU4 zZ_zNc*Qkhb@v>*-5rrr9ftXQE{PT2drG}_MnTk5EM^ddrUBEdzj$7?GL1dYQxGG7`;m@W~f3E;~hDB;^gnq-ecO zXIsTTc~8qbbZVyAQuhSN#!6zyAOm#3Cip4Q{8$VmMipJ-}nV~|CdhpjY z6h$$53pg;1uuU^Y=UP(YQs=3QX>B8MJx2E#Xv=#PPk_7DqZ!}?3LXXLWyPh#r>!|@ z>GI+#Y$Sqa>Mc#Sosc4a`9k;U1^Gpf^3yGuiRo$bxn1(e`;V^_tM&g2@{U}MI^~}xY9v}x;?;)(M`MgyQc}{g zWkYL;;X3ximhzBO;Uh|10(}B0Xe`2K;rt34!)LEAztnxgA#W#uIW5I%o@dE4!(8A3 zG^b|Cv!V!_$jD5OPm$v;A)K2fv6L{`B9C2JW8=&${nsYkQ7#Ur=sAr`OB!l!Ze`Mz zml)5a>g{;bm7O&Wf6%YnUKWe2dfTSAz}*qPtHC!&8M!OovFgL@fn7AM4RdJa@bfBe z{fMiH#}XGcS8-ciKJhDc`NZ!N=i%P~U^_!xNPL^PRIy#F-<1YH2~QAL6Xy^YwNMS_ z5zDt}!Bt2cf^XNtR!r<5K0)l&QpJZnq58{{4BzJ>$WclOiNN^Bn~+wj!e!zJ;uhUh z+(kSP7*|k}-_76!G~)91sy=^Gk>Mugd4=wFml{hbf*SjkK3sp zs>2Z2N(ds>1}TmLmgmSQ z#h(*9#wfmB7bhOpOXUZSRq<)W4&rUZ)^T$Fns$X8h2vC(--xyGiig0+aR3p-0vHcb zdAN$-A`Y3LctjtX5b;J}S%0F6pCqc(EiZ2j*O;+p= zD~I)~rzp-LwoX^PiP%NFpE$xkLnTy$17{RFOR;x<)nMKn#ZKbtc*Wvrl^>C)xC|Hv zP((a9Sj9uKRXm#5D^KwvVtX|?wt)j1$lItAyay=udR6fzU>so3HpPeQc)Q{|#IZXR z*Bhwvi}on)LY%i>@dU+ot?;mNWPw8t;I!guVi&OuPDYHEm#g?r_+uop*LlUJA&QHL zdl75+S5?@e>i8SQ+a=p!#^m^e5{igh4pj|CTvHW>6BiM0sEdEA;>E;XF2z@ft;7$B zT{l#G;xK66j!S+^IV#Wr99yZlM~Gst9~9e(18*zdORW8<*hTF6NwNQMRlocX#a%}z zw%$`bGZe0W?69y#Im{zz0Dmf8HcD~fUy6^xUxG2;Q4dZ8*u0-nTnewt;kJ!9uYuzG z#MO-yZwyoU<#-~&b|Va~e{9gzR5@P?%S|8nC*065-w~tOzLgw>_=*Z_H;GHlQ~~aISiu#iIGb3TtN0V*zy!r^ za7V=a5PY)&wsFL<^A&$b>`0bu*SbZj3WcdE!A4waRa`>s!k0E+t0mU56px;w>IW`V z{33D4BE{#4W0xonn5yy}&mzO|e}No@_lCM%Mh`sPt z6WFd1YtJiwYMQDaLF_kOart_PL;u=taA1Sg=zz^JLvbv=#RA*ZXvL+(-C}3}FROSh zv4i;hOci$#zcWj5$Oe_aY_{SeVC-LuG*g34mC!X-aUSslVr`>}pN~^<2l4JXinUED zz5*UnI3Wk|4EQ*TtmUhCkcHy;O!lvNBq)cMPzl$F3yD9Nr{X~cD*ggI<*>n0;`k)R z5wECtK5^cwioYN(I;6PEd}^tiMNJgXU9Q**-+zN` z!3xC=Z^c`9+Y{>) z7qw9wN9+}-c<}Qo9@tKC%UrnsGzJ~iUO6(rfgDa;K|G(h?+YrvmN#&sc*HJ-S_Z@c! zvIksMQ&fKI^@_`gH?D_ciV5Z9I7|sH;)Z!DUPC;Z*kh_{Fqznk_99M|f5>Iie3WE4~;zHsh#6`p*8&!TO@owTWV$V$~{spo9Avr3@@p`^W zxJkT3D6S^v1&SXM|4!@~rB-D9D=O|wtZ!BvNW6?VNU>e>-=Y#mP{JkRFyeEssrYo_ z!CMu_5^o|-BzAjU#jV6c3Kg#hK^GW4&VB!?FzwB4Z+{E3r^+Yuin?n*r6ZIvHFe1SNY z`1zeGoZj^gKlv471>j+ZFm1K>6u@V$We3bBj0;Vu>TqIU{{h=YhH5ZAqnuzrWJDJ~~&{l4Pr z8FqC{&XYrnR-9X;61<2<>{T2{e4jX&xOkt6hY`;&Ry>`!-G0Rm;-C){yX@rnfgCl& zZyZnwo-t}b(+$Oe#6Bg8Lx|52#}cpqP{kd@nFkdY5kL8nVtXk$d_GnVZKmq*K5-E7 z#!n#L2)^0S`yU61>)s=|OpQ-^A5cel8BaW5~&;M$2EGShK;FWzjpf-mUdl8=}4kX@p zM8$)M7aUa_Qdj?&;t1li#IeMyjzioIOKK&@z!R#1gZM|{JmSKWDqct&S*Ey%c=l(C zONqTdS6ohflGsJO_!QiK@R(MUqt9toL4*H00k=EEUc_&mQSm_HHD4$WA|CLi;t=At z<%%PS`B}xW#QA6O`QKR6@IPz7Z6_okzeIeHxb(J)pC&#=e2ut_xSIGhF+Zo;D<}3Q zK2O~F9Nd4f!wPZ)Q^Hl^vBWOoS;XHHTZ#D{HNcg`fyDX5s?b;o39HoR>;*-S9&Z`;L5_ckQ{j=&In0OHJIN}82*~IgS zZN#aP;rOp2M>-`4;+4d^h@T_=ggBSDoOnI)cf?NOyTpRH-q&hESCOHAEr1+9Qo>Wj zz45;v2-TD7I_EYgNJ?O2{Gplz1)iIpSA{ z|BLuGaSd^yuH!kaHN2!|IKfSEAaNmaU&(O(|4NQ9N*HNU6=o9W6I+QZiB}O1bXWNW z#4ix%Cp7=g-XHQkZ>#~|r4sl1~&xwPH1L~>#@x-ab zvBagsnZR~@`ZcSsDy*f1xx`zD4-gj-dp1z{$B3heD~R_H-yv>*{~8jvD{2B8h`q1S z^S_83T`A!(@nGU}#1X{b5zi(5oj8Zs)JS#s0 z2k{}rc5Nd$yjrLZb`mco{+M{ApNf|guOt4R_zU71;^r+?e#2{OLeq%@iFZkcOiH85@m|Ve+f!laFH2MFI;q&(vIXvk< zBJ=uIB@`0RB(5fYk+_EV6tPDCaKOW*>IW0!{j~{$Hzm%3fvqVr3;_P!foAOApkLDa z!Iyx#x61RHLDU2O1U6mWmtRAM?7DxJaGN+Jz&Mhx?`Sw60rKWSNraLFWdW3AC@D}< zp`<~vK}m;_0VNYk7L;p=^R8pcFuP1P{!BYX{5|qnOCP0}8>a%Ws}jrK(z}d~1`mE*!A}hMkjLjTK5Jv3;4>5-llZ8c z31t?P2T*21F++)k5(lM$E`EAT@81XyJRU1N9&JRg9eOApf{^);Vug|o<$1Auhu*=S z5AroAJD?nZastX(DCeMj1?4;x{4Mkd!~A5ZuM!JGD6@LQm? z2Hg|NQ&9d7$E(f%KVbix17-y5(z|qS8q_Z+xJSQUeVz{Lv7q1l0f~#m)N(!BI8&~B Ml^A!{x~*XU3pnV$e*gdg diff --git a/3D/Utils/mpi_ctvlib.cpp b/3D/Utils/mpi_ctvlib.cpp new file mode 100644 index 0000000..59d6fe4 --- /dev/null +++ b/3D/Utils/mpi_ctvlib.cpp @@ -0,0 +1,574 @@ +// +// tlib.cpp +// TV +// +// Created by Hovden Group on 5/6/19. +// Copyright © 2019 Jonathan Schwartz. All rights reserved. +// + +#include "mpi_ctvlib.hpp" +#include +#include +#include +#include +#include +#include +#include +#include + +#define PI 3.14159265359 + +using namespace Eigen; +using namespace std; +namespace py = pybind11; + +typedef Eigen::Matrix Mat; +typedef Eigen::SparseMatrix SpMat; + +ctvlib::ctvlib(int Ns, int Nray, int Nproj) +{ + //Intialize all the Member variables. + Nslice = Ns; + Ny = Nray; + Nz = Nray; + Nrow = Nray*Nproj; + Ncol = Ny*Nz; + A.resize(Nrow,Ncol); + innerProduct.resize(Nrow); + b.resize(Ny, Nrow); + g.resize(Ny, Nrow); + + MPI_Init(NULL,NULL); + MPI_Comm_size(MPI_COMM_WORLD, &nproc); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + //Calculate the number of slices for each rank. + Nslice_loc = Nslice/nproc; + + //Initialize all the 3D-matrices. + if (rank == 0){ + + //Final Reconstruction. + recon = new Mat[Nslice_loc]; + + // Temporary copy for measuring changes in TV and ART. + temp_recon = new Mat[Nslice_loc]; + + // Temporary copy for measuring 3D TV - Derivative. + tv_recon = new Mat[Nslice_loc]; + + // Original Volume for Simulation Studies. + original_volume = new Mat[Nslice_loc]; + + // Initialize the 3D Matrices as zeros. + #pragma omp parallel for + for (int i=0; i < Nslice_loc; i++) + { + recon[i] = Mat::Zero(Ny, Nz); + temp_recon[i] = Mat::Zero(Ny, Nz); + tv_recon[i] = Mat::Zero(Ny,Nz); + } + } + + //Pass the data to the ranks. + int total_elements = Nslice_loc*Ny*Nz; + MPI_Bcast(&recon, total_elements, MPI_FLOAT, 0, MPI_COMM_WORLD); + MPI_Bcast(&temp_recon, total_elements, MPI_FLOAT, 0, MPI_COMM_WORLD); + MPI_Bcast(&tv_recon, total_elements, MPI_FLOAT, 0, MPI_COMM_WORLD); + +} + +//Import tilt series (projections) from Python. +void ctvlib::setTiltSeries(Mat in) +{ + b = in; +} + +// Import the original volume from python. +void ctvlib::setOriginalVolume(Mat in, int slice) +{ + original_volume[slice] = in; +} + +// Create projections from Volume (for simulation studies) +void ctvlib::create_projections() +{ + #pragma omp parallel for + for (int s = 0; s < Nslice_loc; s++) + { + Mat& mat_slice = original_volume[s]; + mat_slice.resize(mat_slice.size(),1); + VectorXf vec_recon = mat_slice; + for (int i=0; i < Nrow; i++) + { + b(s,i) = A.row(i).dot(vec_recon); + } + mat_slice.resize(Ny,Nz); + } +} + +// Add poisson noise to projections. +void ctvlib::poissonNoise(int Nc) +{ + Mat temp_b = b; + float mean = b.mean(); + float N = b.sum(); + b = b / ( b.sum() ) * Nc * b.size(); + std::default_random_engine generator; + for(int i=0; i < b.size(); i++) + { + std::poisson_distribution distribution(b(i)); + b(i) = distribution(generator); + + } + b = b / ( Nc * b.size() ) * N; + temp_b.array() -= b.array(); + float std = sqrt( ( temp_b.array() - temp_b.mean() ).square().sum() / (temp_b.size() - 1)); +} + +// ART Reconstruction. +void ctvlib::ART(double beta, int dyn_ind) +{ +// + //No dynamic reconstruction, assume fully sampled batch. + if (dyn_ind == -1) { dyn_ind = Nrow; } + //Calculate how many projections were sampled. + else { dyn_ind *= Ny; } + + #pragma omp parallel for + for (int s=0; s < Nslice_loc; s++) + { + Mat& mat_slice = recon[s]; + mat_slice.resize(mat_slice.size(),1); + VectorXf vec_recon = mat_slice; + float a; + for(int j=0; j < dyn_ind; j++) + { + a = ( b(s,j) - A.row(j).dot(vec_recon) ) / innerProduct(j); + vec_recon += A.row(j).transpose() * a * beta; + } + mat_slice = vec_recon; + mat_slice.resize(Ny, Nz); + MPI_Recv(&recon[s], Ny*Ny, MPI_FLOAT, dest, MPI_ANY_TAG, MPI_COMM_WORLD); + } + + MPI_Status status; + if (rank == 0 ) { + right_slice = recon[Nslice_loc]; + MPI_Send(right_slice, Ny*Nz, 1, MPI_FLOAT, rank, MPI_ANY_TAG, MPI_COMM_WORLD); + MPI_Recv(right_slice, Ny*Nz, 1, MPI_FLOAT, rank+1, MPI_ANY_TAG, MPI_COMM_WORLD, &status); + } else if { + right_slice = recon[Nslice_loc]; + MPI_Send(right_slice, Ny*Nz, 1, MPI_FLOAT, rank+1, MPI_ANY_TAG, MPI_COMM_WORLD); + MPI_Recv(right_slice, Ny*Nz, 1, MPI_FLOAT, rank, MPI_ANY_TAG, MPI_COMM_WORLD, &status); + + left_slice = recon[0]; + MPI_Send(left_slice, Ny*Nz, 1, MPI_FLOAT, rank-1, MPI_ANY_TAG, MPI_COMM_WORLD); + MPI_Recv(left_slice, Ny*Nz, 1, MPI_FLOAT, rank, MPI_ANY_TAG, MPI_COMM_WORLD, &status); + } else { + left_slice = recon[0]; + MPI_Send(left_slice, Ny*Nz, 1, MPI_FLOAT, rank-1, MPI_ANY_TAG, MPI_COMM_WORLD); + MPI_Recv(left_slice, Ny*Nz, 1, MPI_FLOAT, rank, MPI_ANY_TAG, MPI_COMM_WORLD, &status); + } +} + +// Stochastic ART Reconstruction. +void ctvlib::sART(double beta, int dyn_ind) +{ + //No dynamic reconstruction, assume fully sampled batch. + if (dyn_ind == -1) { dyn_ind = Nrow; } + //Calculate how many projections were sampled. + else { dyn_ind *= Ny; } + + // Create a random permutation of indices from [0,dyn_ind]. + std::vector r_ind = rand_perm(dyn_ind); + + #pragma omp parallel for + for (int s=0; s < Nslice_loc; s++) + { + Mat& mat_slice = recon[s]; + mat_slice.resize(mat_slice.size(),1); + VectorXf vec_recon = mat_slice; + float a; + for(int j=0; j < dyn_ind; j++) + { + j = r_ind[j]; + a = ( b(s,j) - A.row(j).dot(vec_recon) ) / innerProduct(j); + vec_recon += A.row(j).transpose() * a * beta; + } + mat_slice = vec_recon; + mat_slice.resize(Ny, Nz); + } + + MPI_Status status; + if (rank == 0 ) { + right_slice = recon[Nslice_loc]; + MPI_Send(right_slice, Ny*Nz, 1, MPI_FLOAT, rank, MPI_ANY_TAG, MPI_COMM_WORLD); + MPI_Recv(right_slice, Ny*Nz, 1, MPI_FLOAT, rank+1, MPI_ANY_TAG, MPI_COMM_WORLD, &status); + } else if { + right_slice = recon[Nslice_loc]; + MPI_Send(right_slice, Ny*Nz, 1, MPI_FLOAT, rank+1, MPI_ANY_TAG, MPI_COMM_WORLD); + MPI_Recv(right_slice, Ny*Nz, 1, MPI_FLOAT, rank, MPI_ANY_TAG, MPI_COMM_WORLD, &status); + + left_slice = recon[0]; + MPI_Send(left_slice, Ny*Nz, 1, MPI_FLOAT, rank-1, MPI_ANY_TAG, MPI_COMM_WORLD); + MPI_Recv(left_slice, Ny*Nz, 1, MPI_FLOAT, rank, MPI_ANY_TAG, MPI_COMM_WORLD, &status); + } else { + left_slice = recon[0]; + MPI_Send(left_slice, Ny*Nz, 1, MPI_FLOAT, rank-1, MPI_ANY_TAG, MPI_COMM_WORLD); + MPI_Recv(left_slice, Ny*Nz, 1, MPI_FLOAT, rank, MPI_ANY_TAG, MPI_COMM_WORLD, &status); + } +} + +std::vector ctvlib::rand_perm(int n) +{ + std::vector a(n); + for (int i=0; i < n; i++) + { + a[i] = i; + } + std::random_device rd; + std::mt19937 g(rd()); + std::shuffle(a.begin(), a.end(), g); + return a; +} + +// SIRT Reconstruction. +void ctvlib::SIRT(double beta, int dyn_ind) +{ + //No dynamic reconstruction, assume fully sampled batch. + if (dyn_ind == -1) { dyn_ind = Nrow; } + //Calculate how many projections were sampled. + else { dyn_ind *= Ny; } + + #pragma omp parallel for + for (int s=0; s < Nslice_loc; s++) + { + Mat& mat_slice = recon[s]; + mat_slice.resize(mat_slice.size(),1); + VectorXf vec_recon = mat_slice; + vec_recon += A.transpose() * ( b.row(s).transpose() - A * vec_recon ) * (1/L); + mat_slice = vec_recon; + mat_slice.resize(Ny, Nz); + } + + MPI_Status status; + if (rank == 0 ) { + right_slice = recon[Nslice_loc]; + MPI_Send(right_slice, Ny*Nz, 1, MPI_FLOAT, rank, MPI_ANY_TAG, MPI_COMM_WORLD); + MPI_Recv(right_slice, Ny*Nz, 1, MPI_FLOAT, rank+1, MPI_ANY_TAG, MPI_COMM_WORLD, &status); + } else if { + MPI_Send(right_slice, Ny*Nz, 1, MPI_FLOAT, rank+1, MPI_ANY_TAG, MPI_COMM_WORLD); + MPI_Recv(right_slice, Ny*Nz, 1, MPI_FLOAT, rank, MPI_ANY_TAG, MPI_COMM_WORLD, &status); + + left_slice = recon[0]; + MPI_Send(left_slice, Ny*Nz, 1, MPI_FLOAT, rank-1, MPI_ANY_TAG, MPI_COMM_WORLD); + MPI_Recv(left_slice, Ny*Nz, 1, MPI_FLOAT, rank, MPI_ANY_TAG, MPI_COMM_WORLD, &status); + } else { + left_slice = recon[0]; + MPI_Send(left_slice, Ny*Nz, 1, MPI_FLOAT, rank-1, MPI_ANY_TAG, MPI_COMM_WORLD); + MPI_Recv(left_slice, Ny*Nz, 1, MPI_FLOAT, rank, MPI_ANY_TAG, MPI_COMM_WORLD, &status); + } +} + +//Calculate Lipshits Gradient (for SIRT). +void ctvlib::lipschits() +{ + VectorXf f(Ncol); + f.setOnes(); + float L = (A.transpose() * (A * f)).maxCoeff(); +} + +// Remove Negative Voxels. +void ctvlib::positivity() +{ + #pragma omp parallel for + for(int i=0; i pyA) +{ + for (int i=0; i ctvlib(m, "mpi_ctvlib"); + ctvlib.def(py::init()); + ctvlib.def("setTiltSeries", &ctvlib::setTiltSeries, "Pass the Projections to C++ Object"); + ctvlib.def("setOriginalVolume", &ctvlib::setOriginalVolume, "Pass the Volume to C++ Object"); + ctvlib.def("create_projections", &ctvlib::create_projections, "Create Projections from Volume"); + ctvlib.def("getRecon", &ctvlib::getRecon, "Return the Reconstruction to Python"); + ctvlib.def("ART", &ctvlib::ART, "ART Reconstruction"); + ctvlib.def("sART", &ctvlib::sART, "Stochastic ART Reconstruction"); + ctvlib.def("SIRT", &ctvlib::SIRT, "SIRT Reconstruction"); + ctvlib.def("rowInnerProduct", &ctvlib::normalization, "Calculate the Row Inner Product for Measurement Matrix"); + ctvlib.def("positivity", &ctvlib::positivity, "Remove Negative Elements"); + ctvlib.def("forwardProjection", &ctvlib::forwardProjection, "Forward Projection"); + ctvlib.def("load_A", &ctvlib::loadA, "Load Measurement Matrix Created By Python"); + ctvlib.def("copy_recon", &ctvlib::copy_recon, "Copy the reconstruction"); + ctvlib.def("matrix_2norm", &ctvlib::matrix_2norm, "Calculate L2-Norm of Reconstruction"); + ctvlib.def("vector_2norm", &ctvlib::vector_2norm, "Calculate L2-Norm of Projection (aka Vectors)"); + ctvlib.def("dyn_vector_2norm", &ctvlib::dyn_vector_2norm, "Calculate L2-Norm of Partially Sampled Projections (aka Vectors)"); + ctvlib.def("rmse", &ctvlib::rmse, "Calculate reconstruction's RMSE"); + ctvlib.def("tv", &ctvlib::tv_3D, "Measure 3D TV"); + ctvlib.def("original_tv", &ctvlib::original_tv_3D, "Measure original TV"); + ctvlib.def("tv_gd", &ctvlib::tv_gd_3D, "3D TV Gradient Descent"); + ctvlib.def("get_projections", &ctvlib::get_projections, "Return the projection matrix to python"); + ctvlib.def("poissonNoise", &ctvlib::poissonNoise, "Add Poisson Noise to Projections"); + ctvlib.def("lip", &ctvlib::lipschits, "Add Poisson Noise to Projections"); + ctvlib.def("restart_recon", &ctvlib::restart_recon, "Set all the Slices Equal to Zero"); +} diff --git a/3D/Utils/mpi_ctvlib.hpp b/3D/Utils/mpi_ctvlib.hpp new file mode 100644 index 0000000..7de998c --- /dev/null +++ b/3D/Utils/mpi_ctvlib.hpp @@ -0,0 +1,85 @@ + +// +// mpi_ctvlib.hpp +// TV +// +// Created by Hovden Group on 5/6/19. +// Copyright © 2019 Jonathan Schwartz. All rights reserved. +// + +#ifndef mpi_ctvlib_hpp +#define mpi_ctvlib_hpp + +#include +#include +#include +#include + +class ctvlib +{ + +typedef Eigen::Matrix Mat; +typedef Eigen::SparseMatrix SpMat; + +public: + + // Member Variables. + Mat *recon, *temp_recon, *tv_recon, *original_volume; + Mat& left_slice, right_slice; + SpMat A; + int Nrow, Ncol, Nslice, Nslice_loc, Ny, Nz, nproc, rank; + Eigen::VectorXf innerProduct; + Mat b, g; + + // Initializes Measurement Matrix. + ctvlib(int Nslice, int Nray, int Nproj); + + // Initialize Experimental Projections. + void setTiltSeries(Mat in); + void setOriginalVolume(Mat in, int slice); + void create_projections(); + void poissonNoise(int SNR); + + // Constructs Measurement Matrix. + void loadA(Eigen::Ref pyA); + void normalization(); + + // 2D Reconstructions + void ART(double beta, int dyn_ind); + void SIRT(double beta, int dyn_ind); + void positivity(); + void lipschits(); + + // Stochastic Reconstruction + void sART(double beta, int dyn_ind); + std::vector rand_perm(int n); + + //Forward Project Reconstruction for Data Tolerance Parameter. + void forwardProjection(int dyn_ind); + + // Acquire local copy of reconstruction. + void copy_recon(); + + // Measure 2-norm of projections and reconstruction. + float matrix_2norm(); + float vector_2norm(); + float dyn_vector_2norm(int dyn_ind); + float rmse(); + + // Total variation + float tv_3D(); + float original_tv_3D(); + void tv_gd_3D(int ng, float dPOCS); + + // Set Slices to Zero. + void restart_recon(); + + // Return reconstruction to python. + Mat getRecon(int i); + + // Return projections to python. + Mat get_projections(); + +}; + +#endif /* mpi_ctvlib_hpp */ diff --git a/3D/exp_ASD_tv.py b/3D/exp_ASD_tv.py index 54dc1f7..cfab090 100644 --- a/3D/exp_ASD_tv.py +++ b/3D/exp_ASD_tv.py @@ -11,8 +11,10 @@ import time ######################################## -vol_size = '256_' -file_name = 'Co2P_tiltser.tif' +# vol_size = '256_' +# file_name = '180_tiltser.tif' +vol_size = '' +file_name = 'FePt_tiltser.npy' # Number of Iterations (Main Loop) Niter = 100 @@ -27,7 +29,7 @@ beta_red = 0.99 # Data Tolerance Parameter -eps = 0.43 +eps = 0.002 # Reduction Criteria r_max = 0.95 diff --git a/3D/sim_tv_mpi.py b/3D/sim_tv_mpi.py new file mode 100644 index 0000000..88eec28 --- /dev/null +++ b/3D/sim_tv_mpi.py @@ -0,0 +1,157 @@ +# General 3D - ASD/TV Reconstruction with Positivity Constraint. +# Intended for simulated datasets to measure RMSE and Volume's Original TV. +# and to reconstruct large volume sizes (>1000^3) with Distributed Memory (OpenMPI) + +import sys, os +sys.path.append('../Utils') +from pytvlib import parallelRay, timer, load_data +import plot_results as pr +import numpy as np +import mpi_ctvlib +import time +######################################## + +vol_size = '256_' +file_name = 'au_sto_tiltser.npy' + +# Number of Iterations (Main Loop) +Niter = 300 + +# Number of Iterations (TV Loop) +ng = 10 + +# Parameter in ART Reconstruction. +beta = 0.25 + +# ART Reduction. +beta_red = 0.985 + +# Data Tolerance Parameter +eps = 0.019 + +# Reduction Criteria +r_max = 0.95 +alpha_red = 0.95 +alpha = 0.2 + +SNR = 100 + +#Outcomes: +noise = True # Add noise to the reconstruction. +show_live_plot = 0 +save_recon = 0 # Save final Reconstruction. +########################################## + +#Read Image. +(file_name, original_volume) = load_data(vol_size,file_name) +file_name = 'au_sto' +(Nslice, Nray, _) = original_volume.shape + +# Generate Tilt Angles. +tiltAngles = np.load('Tilt_Series/'+ file_name +'_tiltAngles.npy') +Nproj = tiltAngles.shape[0] + +# Initialize C++ Object.. +tomo_obj = ctvlib.ctvlib(Nslice, Nray, Nproj) + +# Generate measurement matrix +A = parallelRay(Nray, tiltAngles) +tomo_obj.load_A(A) +A = None +tomo_obj.rowInnerProduct() + +print('Measurement Matrix is Constructed!') + +# If creating simulation with noise, set background value to 1. +if noise: + original_volume[original_volume == 0] = 1 + +# Load Volume and Collect Projections. +for s in range(Nslice): + tomo_obj.setOriginalVolume(original_volume[s,:,:],s) +tomo_obj.create_projections() + +# Apply poisson noise to volume. +if noise: + tomo_obj.poissonNoise(SNR) + +#Measure Volume's Original TV +tv0 = tomo_obj.original_tv() + +gif = np.zeros([Nray, Nray, Niter], dtype=np.float32) +gif2 = np.zeros([Nray, Nray, Niter], dtype=np.float32) + +dd_vec = np.zeros(Niter) +tv_vec = np.zeros(Niter) +rmse_vec = np.zeros(Niter) +time_vec = np.zeros(Niter) + +counter = 1 + +t0 = time.time() + +#Main Loop +for i in range(Niter): + + if ( i % 25 ==0): + print('Iteration No.: ' + str(i+1) +'/'+str(Niter)) + + tomo_obj.copy_recon() + + #ART Reconstruction. + tomo_obj.sART(beta, -1) + + #Positivity constraint + tomo_obj.positivity() + + #ART-Beta Reduction + beta *= beta_red + + #Forward Projection. + tomo_obj.forwardProjection(-1) + + #Measure Magnitude for TV - GD. + if (i == 0): + dPOCS = tomo_obj.matrix_2norm() * alpha + dp = dPOCS / alpha + else: # Measure change from ART. + dp = tomo_obj.matrix_2norm() + + # Measure difference between exp/sim projections. + dd_vec[i] = tomo_obj.vector_2norm() + + #Measure TV. + tv_vec[i] = tomo_obj.tv() + + #Measure RMSE. + rmse_vec[i] = tomo_obj.rmse() + + tomo_obj.copy_recon() + + #TV Minimization. + tomo_obj.tv_gd(ng, dPOCS) + dg = tomo_obj.matrix_2norm() + + if(dg > dp * r_max and dd_vec[i] > eps): + dPOCS *= alpha_red + + if (i+1)% 25 == 0: + timer(t0, counter, Niter) + if show_live_plot: + pr.sim_ASD_live_plot(dd_vec, eps, tv_vec, tv0, rmse_vec, i) + counter += 1 + time_vec[i] = time.time() - t0 + +#Save all the results to single matrix. +results = np.array([dd_vec, eps, tv_vec, tv0, rmse_vec, time_vec]) +os.makedirs('Results/'+ file_name +'_ASD/', exist_ok=True) +np.save('Results/' + file_name + '_ASD/results5.npy', results) + +#Get and save the final reconstruction. +if save_recon: + recon = np.zeros([Nslice, Nray, Nray], dtype=np.float32, order='F') + + for s in range(Nslice): + recon[s,:,:] = tomo_obj.getRecon(s) + + np.save('Results/TV_'+ file_name + '_recon.npy', recon)