From d54172290f629d98c8a330ad358d5932bdba8dbc Mon Sep 17 00:00:00 2001 From: Jonas Adler Date: Tue, 19 Apr 2016 09:30:20 +0200 Subject: Fix missing backslash in readme.md install --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index b34cf38..fb1f641 100644 --- a/README.md +++ b/README.md @@ -30,7 +30,7 @@ cd build/linux ./autogen.sh # when building a git version ./configure --with-cuda=/usr/local/cuda \ --with-matlab=/usr/local/MATLAB/R2012a \ - --with-python + --with-python \ --prefix=/usr/local/astra make make install -- cgit v1.2.3 From c58929ca7d0201335b5036773b26f82361793f25 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Tue, 19 Apr 2016 16:40:14 +0200 Subject: Use CXX for Python compilation as well --- build/linux/Makefile.in | 4 ++-- python/astra/extrautils.pyx | 2 ++ 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in index a199bf6..2214734 100644 --- a/build/linux/Makefile.in +++ b/build/linux/Makefile.in @@ -281,13 +281,13 @@ ifeq ($(python),yes) py: libastra.la $(MKDIR) python/build $(MKDIR) python/finalbuild - cd $(srcdir)/../../python; CPPFLAGS="${PYCPPFLAGS}" LDFLAGS="${PYLDFLAGS}" $(PYTHON) builder.py build --build-base=$(abs_top_builddir)/python/build install \ + cd $(srcdir)/../../python; CXX="${CXX}" CC="${CXX}" CPPFLAGS="${PYCPPFLAGS}" LDFLAGS="${PYLDFLAGS}" $(PYTHON) builder.py build --build-base=$(abs_top_builddir)/python/build install \ --install-base=$(abs_top_builddir)/python/finalbuild --install-headers=$(abs_top_builddir)/python/finalbuild --install-purelib=$(abs_top_builddir)/python/finalbuild \ --install-platlib=$(abs_top_builddir)/python/finalbuild --install-scripts=$(abs_top_builddir)/python/finalbuild --install-data=$(abs_top_builddir)/python/finalbuild python-root-install: libastra.la $(MKDIR) python/build - cd $(srcdir)/../../python; CPPFLAGS="${PYCPPFLAGS}" LDFLAGS="${PYLDFLAGS}" $(PYTHON) builder.py build --build-base=$(abs_top_builddir)/python/build install + cd $(srcdir)/../../python; CXX="${CXX}" CC="${CXX}" CPPFLAGS="${PYCPPFLAGS}" LDFLAGS="${PYLDFLAGS}" $(PYTHON) builder.py build --build-base=$(abs_top_builddir)/python/build install endif diff --git a/python/astra/extrautils.pyx b/python/astra/extrautils.pyx index 5bc315f..2c7771e 100644 --- a/python/astra/extrautils.pyx +++ b/python/astra/extrautils.pyx @@ -22,6 +22,8 @@ # along with the ASTRA Toolbox. If not, see . # # ----------------------------------------------------------------------- +# distutils: language = c++ + def clipCircle(img): cdef int i,j -- cgit v1.2.3 From 93eacadc5b94f8e1241dc6836be27d16881a2826 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Tue, 19 Apr 2016 17:12:58 +0200 Subject: Add comment about setting CC to CXX --- build/linux/Makefile.in | 2 ++ 1 file changed, 2 insertions(+) diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in index 2214734..f10f482 100644 --- a/build/linux/Makefile.in +++ b/build/linux/Makefile.in @@ -281,12 +281,14 @@ ifeq ($(python),yes) py: libastra.la $(MKDIR) python/build $(MKDIR) python/finalbuild + # Note: setting CC to CXX is intentional. Python uses CC for compilation even if input is C++. cd $(srcdir)/../../python; CXX="${CXX}" CC="${CXX}" CPPFLAGS="${PYCPPFLAGS}" LDFLAGS="${PYLDFLAGS}" $(PYTHON) builder.py build --build-base=$(abs_top_builddir)/python/build install \ --install-base=$(abs_top_builddir)/python/finalbuild --install-headers=$(abs_top_builddir)/python/finalbuild --install-purelib=$(abs_top_builddir)/python/finalbuild \ --install-platlib=$(abs_top_builddir)/python/finalbuild --install-scripts=$(abs_top_builddir)/python/finalbuild --install-data=$(abs_top_builddir)/python/finalbuild python-root-install: libastra.la $(MKDIR) python/build + # Note: setting CC to CXX is intentional. Python uses CC for compilation even if input is C++. cd $(srcdir)/../../python; CXX="${CXX}" CC="${CXX}" CPPFLAGS="${PYCPPFLAGS}" LDFLAGS="${PYLDFLAGS}" $(PYTHON) builder.py build --build-base=$(abs_top_builddir)/python/build install endif -- cgit v1.2.3 From 7f8b7b37121e06cec76b2b1cab1b56920c2c4ef4 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 22 Apr 2016 17:15:03 +0200 Subject: Fix build --- include/astra/Fourier.h | 2 +- src/FilteredBackProjectionAlgorithm.cpp | 2 +- src/Fourier.cpp | 206 ++++++++++++-------------------- 3 files changed, 78 insertions(+), 132 deletions(-) diff --git a/include/astra/Fourier.h b/include/astra/Fourier.h index ff26f82..68f9f38 100644 --- a/include/astra/Fourier.h +++ b/include/astra/Fourier.h @@ -76,7 +76,7 @@ namespace astra { } . */ -void cdft(int n, int isgn, float32 *a, int *ip, float32 *w); +_AstraExport void cdft(int n, int isgn, float32 *a, int *ip, float32 *w); } diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp index 90efd52..70462f7 100644 --- a/src/FilteredBackProjectionAlgorithm.cpp +++ b/src/FilteredBackProjectionAlgorithm.cpp @@ -275,7 +275,7 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D float32* pf = new float32[2 * iAngleCount * zpDetector]; - int *ip = new int[int(2+sqrt(zpDetector)+1)]; + int *ip = new int[int(2+sqrt((float)zpDetector)+1)]; ip[0]=0; float32 *w = new float32[zpDetector/2]; diff --git a/src/Fourier.cpp b/src/Fourier.cpp index 5ca22e6..c33f7bd 100644 --- a/src/Fourier.cpp +++ b/src/Fourier.cpp @@ -27,6 +27,7 @@ $Id$ */ #include "astra/Fourier.h" +#include namespace astra { @@ -320,11 +321,45 @@ Appendix : */ -void cdft(int n, int isgn, float32 *a, int *ip, float32 *w) +static int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w); +static void bitrv208(float32 *a); +static void bitrv208neg(float32 *a); +static void bitrv216(float32 *a); +static void bitrv216neg(float32 *a); +static void bitrv2conj(int n, int *ip, float32 *a); +static void bitrv2(int n, int *ip, float32 *a); +static void cftb040(float32 *a); +static void cftb1st(int n, float32 *a, float32 *w); +static void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w); +static void cftf040(float32 *a); +static void cftf081(float32 *a, float32 *w); +static void cftf082(float32 *a, float32 *w); +static void cftf161(float32 *a, float32 *w); +static void cftf162(float32 *a, float32 *w); +static void cftf1st(int n, float32 *a, float32 *w); +static void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w); +static void cftfx41(int n, float32 *a, int nw, float32 *w); +static void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w); +static void cftmdl1(int n, float32 *a, float32 *w); +static void cftmdl2(int n, float32 *a, float32 *w); +static void *cftrec1_th(void *p); +static void *cftrec2_th(void *p); +static void cftrec4(int n, float32 *a, int nw, float32 *w); +static void cftx020(float32 *a); +static void dctsub(int n, float32 *a, int nc, float32 *c); +static void dstsub(int n, float32 *a, int nc, float32 *c); +static void makect(int nc, int *ip, float32 *c); +static void makeipt(int nw, int *ip); +static void makewt(int nw, int *ip, float32 *w); +static void rftbsub(int n, float32 *a, int nc, float32 *c); +static void rftfsub(int n, float32 *a, int nc, float32 *c); +#ifdef USE_CDFT_THREADS +static void cftrec4_th(int n, float32 *a, int nw, float32 *w); +#endif /* USE_CDFT_THREADS */ + + +_AstraExport void cdft(int n, int isgn, float32 *a, int *ip, float32 *w) { - void makewt(int nw, int *ip, float32 *w); - void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w); - void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w); int nw; nw = ip[0]; @@ -340,14 +375,8 @@ void cdft(int n, int isgn, float32 *a, int *ip, float32 *w) } -void rdft(int n, int isgn, float32 *a, int *ip, float32 *w) +_AstraExport void rdft(int n, int isgn, float32 *a, int *ip, float32 *w) { - void makewt(int nw, int *ip, float32 *w); - void makect(int nc, int *ip, float32 *c); - void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w); - void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w); - void rftfsub(int n, float32 *a, int nc, float32 *c); - void rftbsub(int n, float32 *a, int nc, float32 *c); int nw, nc; float32 xi; @@ -384,15 +413,8 @@ void rdft(int n, int isgn, float32 *a, int *ip, float32 *w) } -void ddct(int n, int isgn, float32 *a, int *ip, float32 *w) +_AstraExport void ddct(int n, int isgn, float32 *a, int *ip, float32 *w) { - void makewt(int nw, int *ip, float32 *w); - void makect(int nc, int *ip, float32 *c); - void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w); - void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w); - void rftfsub(int n, float32 *a, int nc, float32 *c); - void rftbsub(int n, float32 *a, int nc, float32 *c); - void dctsub(int n, float32 *a, int nc, float32 *c); int j, nw, nc; float32 xr; @@ -440,15 +462,8 @@ void ddct(int n, int isgn, float32 *a, int *ip, float32 *w) } -void ddst(int n, int isgn, float32 *a, int *ip, float32 *w) +_AstraExport void ddst(int n, int isgn, float32 *a, int *ip, float32 *w) { - void makewt(int nw, int *ip, float32 *w); - void makect(int nc, int *ip, float32 *c); - void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w); - void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w); - void rftfsub(int n, float32 *a, int nc, float32 *c); - void rftbsub(int n, float32 *a, int nc, float32 *c); - void dstsub(int n, float32 *a, int nc, float32 *c); int j, nw, nc; float32 xr; @@ -496,13 +511,8 @@ void ddst(int n, int isgn, float32 *a, int *ip, float32 *w) } -void dfct(int n, float32 *a, float32 *t, int *ip, float32 *w) +_AstraExport void dfct(int n, float32 *a, float32 *t, int *ip, float32 *w) { - void makewt(int nw, int *ip, float32 *w); - void makect(int nc, int *ip, float32 *c); - void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w); - void rftfsub(int n, float32 *a, int nc, float32 *c); - void dctsub(int n, float32 *a, int nc, float32 *c); int j, k, l, m, mh, nw, nc; float32 xr, xi, yr, yi; @@ -589,13 +599,8 @@ void dfct(int n, float32 *a, float32 *t, int *ip, float32 *w) } -void dfst(int n, float32 *a, float32 *t, int *ip, float32 *w) +_AstraExport void dfst(int n, float32 *a, float32 *t, int *ip, float32 *w) { - void makewt(int nw, int *ip, float32 *w); - void makect(int nc, int *ip, float32 *c); - void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w); - void rftfsub(int n, float32 *a, int nc, float32 *c); - void dstsub(int n, float32 *a, int nc, float32 *c); int j, k, l, m, mh, nw, nc; float32 xr, xi, yr, yi; @@ -675,12 +680,8 @@ void dfst(int n, float32 *a, float32 *t, int *ip, float32 *w) /* -------- initializing routines -------- */ - -#include - -void makewt(int nw, int *ip, float32 *w) +static void makewt(int nw, int *ip, float32 *w) { - void makeipt(int nw, int *ip); int j, nwh, nw0, nw1; float32 delta, wn4r, wk1r, wk1i, wk3r, wk3i; @@ -739,7 +740,7 @@ void makewt(int nw, int *ip, float32 *w) } -void makeipt(int nw, int *ip) +static void makeipt(int nw, int *ip) { int j, l, m, m2, p, q; @@ -759,7 +760,7 @@ void makeipt(int nw, int *ip) } -void makect(int nc, int *ip, float32 *c) +static void makect(int nc, int *ip, float32 *c) { int j, nch; float32 delta; @@ -835,23 +836,8 @@ void makect(int nc, int *ip, float32 *c) #endif /* USE_CDFT_WINTHREADS */ -void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w) +static void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w) { - void bitrv2(int n, int *ip, float32 *a); - void bitrv216(float32 *a); - void bitrv208(float32 *a); - void cftf1st(int n, float32 *a, float32 *w); - void cftrec4(int n, float32 *a, int nw, float32 *w); - void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w); - void cftfx41(int n, float32 *a, int nw, float32 *w); - void cftf161(float32 *a, float32 *w); - void cftf081(float32 *a, float32 *w); - void cftf040(float32 *a); - void cftx020(float32 *a); -#ifdef USE_CDFT_THREADS - void cftrec4_th(int n, float32 *a, int nw, float32 *w); -#endif /* USE_CDFT_THREADS */ - if (n > 8) { if (n > 32) { cftf1st(n, a, &w[nw - (n >> 2)]); @@ -883,23 +869,8 @@ void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w) } -void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w) +static void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w) { - void bitrv2conj(int n, int *ip, float32 *a); - void bitrv216neg(float32 *a); - void bitrv208neg(float32 *a); - void cftb1st(int n, float32 *a, float32 *w); - void cftrec4(int n, float32 *a, int nw, float32 *w); - void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w); - void cftfx41(int n, float32 *a, int nw, float32 *w); - void cftf161(float32 *a, float32 *w); - void cftf081(float32 *a, float32 *w); - void cftb040(float32 *a); - void cftx020(float32 *a); -#ifdef USE_CDFT_THREADS - void cftrec4_th(int n, float32 *a, int nw, float32 *w); -#endif /* USE_CDFT_THREADS */ - if (n > 8) { if (n > 32) { cftb1st(n, a, &w[nw - (n >> 2)]); @@ -931,7 +902,7 @@ void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w) } -void bitrv2(int n, int *ip, float32 *a) +static void bitrv2(int n, int *ip, float32 *a) { int j, j1, k, k1, l, m, nh, nm; float32 xr, xi, yr, yi; @@ -1278,7 +1249,7 @@ void bitrv2(int n, int *ip, float32 *a) } -void bitrv2conj(int n, int *ip, float32 *a) +static void bitrv2conj(int n, int *ip, float32 *a) { int j, j1, k, k1, l, m, nh, nm; float32 xr, xi, yr, yi; @@ -1633,7 +1604,7 @@ void bitrv2conj(int n, int *ip, float32 *a) } -void bitrv216(float32 *a) +static void bitrv216(float32 *a) { float32 x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i, @@ -1690,7 +1661,7 @@ void bitrv216(float32 *a) } -void bitrv216neg(float32 *a) +static void bitrv216neg(float32 *a) { float32 x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i, @@ -1760,7 +1731,7 @@ void bitrv216neg(float32 *a) } -void bitrv208(float32 *a) +static void bitrv208(float32 *a) { float32 x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i; @@ -1783,7 +1754,7 @@ void bitrv208(float32 *a) } -void bitrv208neg(float32 *a) +static void bitrv208neg(float32 *a) { float32 x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, x5r, x5i, x6r, x6i, x7r, x7i; @@ -1819,7 +1790,7 @@ void bitrv208neg(float32 *a) } -void cftf1st(int n, float32 *a, float32 *w) +static void cftf1st(int n, float32 *a, float32 *w) { int j, j0, j1, j2, j3, k, m, mh; float32 wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, @@ -2025,7 +1996,7 @@ void cftf1st(int n, float32 *a, float32 *w) } -void cftb1st(int n, float32 *a, float32 *w) +static void cftb1st(int n, float32 *a, float32 *w) { int j, j0, j1, j2, j3, k, m, mh; float32 wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, @@ -2242,10 +2213,8 @@ struct cdft_arg_st { typedef struct cdft_arg_st cdft_arg_t; -void cftrec4_th(int n, float32 *a, int nw, float32 *w) +static void cftrec4_th(int n, float32 *a, int nw, float32 *w) { - void *cftrec1_th(void *p); - void *cftrec2_th(void *p); int i, idiv4, m, nthread; cdft_thread_t th[4]; cdft_arg_t ag[4]; @@ -2276,11 +2245,8 @@ void cftrec4_th(int n, float32 *a, int nw, float32 *w) } -void *cftrec1_th(void *p) +static void *cftrec1_th(void *p) { - int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w); - void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w); - void cftmdl1(int n, float32 *a, float32 *w); int isplt, j, k, m, n, n0, nw; float32 *a, *w; @@ -2305,11 +2271,8 @@ void *cftrec1_th(void *p) } -void *cftrec2_th(void *p) +static void *cftrec2_th(void *p) { - int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w); - void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w); - void cftmdl2(int n, float32 *a, float32 *w); int isplt, j, k, m, n, n0, nw; float32 *a, *w; @@ -2337,11 +2300,8 @@ void *cftrec2_th(void *p) #endif /* USE_CDFT_THREADS */ -void cftrec4(int n, float32 *a, int nw, float32 *w) +static void cftrec4(int n, float32 *a, int nw, float32 *w) { - int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w); - void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w); - void cftmdl1(int n, float32 *a, float32 *w); int isplt, j, k, m; m = n; @@ -2361,8 +2321,6 @@ void cftrec4(int n, float32 *a, int nw, float32 *w) int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w) { - void cftmdl1(int n, float32 *a, float32 *w); - void cftmdl2(int n, float32 *a, float32 *w); int i, isplt, m; if ((k & 3) != 0) { @@ -2394,15 +2352,8 @@ int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w) } -void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w) +static void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w) { - void cftmdl1(int n, float32 *a, float32 *w); - void cftmdl2(int n, float32 *a, float32 *w); - void cftf161(float32 *a, float32 *w); - void cftf162(float32 *a, float32 *w); - void cftf081(float32 *a, float32 *w); - void cftf082(float32 *a, float32 *w); - if (n == 512) { cftmdl1(128, a, &w[nw - 64]); cftf161(a, &w[nw - 8]); @@ -2459,7 +2410,7 @@ void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w) } -void cftmdl1(int n, float32 *a, float32 *w) +static void cftmdl1(int n, float32 *a, float32 *w) { int j, j0, j1, j2, j3, k, m, mh; float32 wn4r, wk1r, wk1i, wk3r, wk3i; @@ -2569,7 +2520,7 @@ void cftmdl1(int n, float32 *a, float32 *w) } -void cftmdl2(int n, float32 *a, float32 *w) +static void cftmdl2(int n, float32 *a, float32 *w) { int j, j0, j1, j2, j3, k, kr, m, mh; float32 wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i; @@ -2703,13 +2654,8 @@ void cftmdl2(int n, float32 *a, float32 *w) } -void cftfx41(int n, float32 *a, int nw, float32 *w) +static void cftfx41(int n, float32 *a, int nw, float32 *w) { - void cftf161(float32 *a, float32 *w); - void cftf162(float32 *a, float32 *w); - void cftf081(float32 *a, float32 *w); - void cftf082(float32 *a, float32 *w); - if (n == 128) { cftf161(a, &w[nw - 8]); cftf162(&a[32], &w[nw - 32]); @@ -2724,7 +2670,7 @@ void cftfx41(int n, float32 *a, int nw, float32 *w) } -void cftf161(float32 *a, float32 *w) +static void cftf161(float32 *a, float32 *w) { float32 wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, @@ -2883,7 +2829,7 @@ void cftf161(float32 *a, float32 *w) } -void cftf162(float32 *a, float32 *w) +static void cftf162(float32 *a, float32 *w) { float32 wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, x0r, x0i, x1r, x1i, x2r, x2i, @@ -3066,7 +3012,7 @@ void cftf162(float32 *a, float32 *w) } -void cftf081(float32 *a, float32 *w) +static void cftf081(float32 *a, float32 *w) { float32 wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, @@ -3128,7 +3074,7 @@ void cftf081(float32 *a, float32 *w) } -void cftf082(float32 *a, float32 *w) +static void cftf082(float32 *a, float32 *w) { float32 wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, @@ -3200,7 +3146,7 @@ void cftf082(float32 *a, float32 *w) } -void cftf040(float32 *a) +static void cftf040(float32 *a) { float32 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; @@ -3223,7 +3169,7 @@ void cftf040(float32 *a) } -void cftb040(float32 *a) +static void cftb040(float32 *a) { float32 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; @@ -3246,7 +3192,7 @@ void cftb040(float32 *a) } -void cftx020(float32 *a) +static void cftx020(float32 *a) { float32 x0r, x0i; @@ -3259,7 +3205,7 @@ void cftx020(float32 *a) } -void rftfsub(int n, float32 *a, int nc, float32 *c) +static void rftfsub(int n, float32 *a, int nc, float32 *c) { int j, k, kk, ks, m; float32 wkr, wki, xr, xi, yr, yi; @@ -3284,7 +3230,7 @@ void rftfsub(int n, float32 *a, int nc, float32 *c) } -void rftbsub(int n, float32 *a, int nc, float32 *c) +static void rftbsub(int n, float32 *a, int nc, float32 *c) { int j, k, kk, ks, m; float32 wkr, wki, xr, xi, yr, yi; @@ -3309,7 +3255,7 @@ void rftbsub(int n, float32 *a, int nc, float32 *c) } -void dctsub(int n, float32 *a, int nc, float32 *c) +static void dctsub(int n, float32 *a, int nc, float32 *c) { int j, k, kk, ks, m; float32 wkr, wki, xr; @@ -3330,7 +3276,7 @@ void dctsub(int n, float32 *a, int nc, float32 *c) } -void dstsub(int n, float32 *a, int nc, float32 *c) +static void dstsub(int n, float32 *a, int nc, float32 *c) { int j, k, kk, ks, m; float32 wkr, wki, xr; -- cgit v1.2.3 From 8b67986464daae799d0171aed70a0d2cd96fd8d1 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 22 Apr 2016 17:39:59 +0200 Subject: Fix build --- include/astra/AstraObjectManager.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/astra/AstraObjectManager.h b/include/astra/AstraObjectManager.h index ad89c2a..9faecbe 100644 --- a/include/astra/AstraObjectManager.h +++ b/include/astra/AstraObjectManager.h @@ -60,7 +60,7 @@ public: }; -class CAstraIndexManager : public Singleton { +class _AstraExport CAstraIndexManager : public Singleton { public: CAstraIndexManager() : m_iLastIndex(0) { } -- cgit v1.2.3 From 1e26f7602b6685c584fd4d857353f390622e3a34 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Mon, 25 Apr 2016 10:47:59 +0200 Subject: Change flatten to ravel in Python code --- python/astra/optomo.py | 4 ++-- samples/python/s009_projection_matrix.py | 2 +- samples/python/s015_fp_bp.py | 6 +++--- samples/python/s017_OpTomo.py | 2 +- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/python/astra/optomo.py b/python/astra/optomo.py index dd10713..5a92998 100644 --- a/python/astra/optomo.py +++ b/python/astra/optomo.py @@ -125,7 +125,7 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): algorithm.delete(fp_id) self.data_mod.delete([vid,sid]) - return s.flatten() + return s.ravel() def rmatvec(self,s): """Implements the transpose operator. @@ -147,7 +147,7 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): algorithm.delete(bp_id) self.data_mod.delete([vid,sid]) - return v.flatten() + return v.ravel() def __mul__(self,v): """Provides easy forward operator by *. diff --git a/samples/python/s009_projection_matrix.py b/samples/python/s009_projection_matrix.py index c4c4557..e20d58c 100644 --- a/samples/python/s009_projection_matrix.py +++ b/samples/python/s009_projection_matrix.py @@ -46,7 +46,7 @@ W = astra.matrix.get(matrix_id) # Manually use this projection matrix to do a projection: import scipy.io P = scipy.io.loadmat('phantom.mat')['phantom256'] -s = W.dot(P.flatten()) +s = W.dot(P.ravel()) s = np.reshape(s, (len(proj_geom['ProjectionAngles']),proj_geom['DetectorCount'])) import pylab diff --git a/samples/python/s015_fp_bp.py b/samples/python/s015_fp_bp.py index fa0bf86..ff0b30a 100644 --- a/samples/python/s015_fp_bp.py +++ b/samples/python/s015_fp_bp.py @@ -46,12 +46,12 @@ class astra_wrap(object): def matvec(self,v): sid, s = astra.create_sino(np.reshape(v,(vol_geom['GridRowCount'],vol_geom['GridColCount'])),self.proj_id) astra.data2d.delete(sid) - return s.flatten() + return s.ravel() def rmatvec(self,v): bid, b = astra.create_backprojection(np.reshape(v,(len(proj_geom['ProjectionAngles']),proj_geom['DetectorCount'],)),self.proj_id) astra.data2d.delete(bid) - return b.flatten() + return b.ravel() vol_geom = astra.create_vol_geom(256, 256) proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False)) @@ -65,7 +65,7 @@ proj_id = astra.create_projector('cuda',proj_geom,vol_geom) sinogram_id, sinogram = astra.create_sino(P, proj_id) # Reshape the sinogram into a vector -b = sinogram.flatten() +b = sinogram.ravel() # Call lsqr with ASTRA FP and BP import scipy.sparse.linalg diff --git a/samples/python/s017_OpTomo.py b/samples/python/s017_OpTomo.py index 967fa64..214e9a7 100644 --- a/samples/python/s017_OpTomo.py +++ b/samples/python/s017_OpTomo.py @@ -50,7 +50,7 @@ pylab.figure(2) pylab.imshow(sinogram) # Run the lsqr linear solver -output = scipy.sparse.linalg.lsqr(W, sinogram.flatten(), iter_lim=150) +output = scipy.sparse.linalg.lsqr(W, sinogram.ravel(), iter_lim=150) rec = output[0].reshape([256, 256]) pylab.figure(3) -- cgit v1.2.3 From a33b2a0e4e0663145fb0d4d67a9145948f90341a Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 25 Apr 2016 11:03:18 +0200 Subject: Reverse output of astra.data3d.dimensions This makes the output directly usable in numpy commands like zeros, reshape. We were not using data3d.dimensions ourselves, as far as I can tell. --- python/astra/data3d_c.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx index 207d9a5..811d1e4 100644 --- a/python/astra/data3d_c.pyx +++ b/python/astra/data3d_c.pyx @@ -264,7 +264,7 @@ def store(i,data): def dimensions(i): cdef CFloat32Data3D * pDataObject = getObject(i) - return (pDataObject.getWidth(),pDataObject.getHeight(),pDataObject.getDepth()) + return (pDataObject.getDepth(),pDataObject.getHeight(),pDataObject.getWidth()) def delete(ids): try: -- cgit v1.2.3 From 63e22423350089a37b188a53164867eedab781ac Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 21 Apr 2016 11:39:51 +0200 Subject: Give OpTomo FP/BP functions with optional out argument This allows more efficient use of allocated arrays. --- python/astra/optomo.py | 96 ++++++++++++++++++++++++++++++++++---------------- 1 file changed, 66 insertions(+), 30 deletions(-) diff --git a/python/astra/optomo.py b/python/astra/optomo.py index 5a92998..dde719e 100644 --- a/python/astra/optomo.py +++ b/python/astra/optomo.py @@ -111,21 +111,7 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): :param v: Volume to forward project. :type v: :class:`numpy.ndarray` """ - v = self.__checkArray(v, self.vshape) - vid = self.data_mod.link('-vol',self.vg,v) - s = np.zeros(self.sshape,dtype=np.float32) - sid = self.data_mod.link('-sino',self.pg,s) - - cfg = creators.astra_dict('FP'+self.appendString) - cfg['ProjectionDataId'] = sid - cfg['VolumeDataId'] = vid - cfg['ProjectorId'] = self.proj_id - fp_id = algorithm.create(cfg) - algorithm.run(fp_id) - - algorithm.delete(fp_id) - self.data_mod.delete([vid,sid]) - return s.ravel() + return self.FP(v, out=None).ravel() def rmatvec(self,s): """Implements the transpose operator. @@ -133,21 +119,7 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): :param s: The projection data. :type s: :class:`numpy.ndarray` """ - s = self.__checkArray(s, self.sshape) - sid = self.data_mod.link('-sino',self.pg,s) - v = np.zeros(self.vshape,dtype=np.float32) - vid = self.data_mod.link('-vol',self.vg,v) - - cfg = creators.astra_dict('BP'+self.appendString) - cfg['ProjectionDataId'] = sid - cfg['ReconstructionDataId'] = vid - cfg['ProjectorId'] = self.proj_id - bp_id = algorithm.create(cfg) - algorithm.run(bp_id) - - algorithm.delete(bp_id) - self.data_mod.delete([vid,sid]) - return v.ravel() + return self.BP(s, out=None).ravel() def __mul__(self,v): """Provides easy forward operator by *. @@ -189,6 +161,70 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): self.data_mod.delete([vid,sid]) return v + def FP(self,v,out=None): + """Perform forward projection. + + Output must have the right 2D/3D shape. Input may also be flattened. + + Output must also be contiguous and float32. This isn't required for the + input, but it is more efficient if it is. + + :param v: Volume to forward project. + :type v: :class:`numpy.ndarray` + :param out: Array to store result in. + :type out: :class:`numpy.ndarray` + """ + + v = self.__checkArray(v, self.vshape) + vid = self.data_mod.link('-vol',self.vg,v) + if out is None: + out = np.zeros(self.sshape,dtype=np.float32) + sid = self.data_mod.link('-sino',self.pg,out) + + cfg = creators.astra_dict('FP'+self.appendString) + cfg['ProjectionDataId'] = sid + cfg['VolumeDataId'] = vid + cfg['ProjectorId'] = self.proj_id + fp_id = algorithm.create(cfg) + algorithm.run(fp_id) + + algorithm.delete(fp_id) + self.data_mod.delete([vid,sid]) + return out + + def BP(self,s,out=None): + """Perform backprojection. + + Output must have the right 2D/3D shape. Input may also be flattened. + + Output must also be contiguous and float32. This isn't required for the + input, but it is more efficient if it is. + + :param : The projection data. + :type s: :class:`numpy.ndarray` + :param out: Array to store result in. + :type out: :class:`numpy.ndarray` + """ + s = self.__checkArray(s, self.sshape) + sid = self.data_mod.link('-sino',self.pg,s) + if out is None: + out = np.zeros(self.vshape,dtype=np.float32) + vid = self.data_mod.link('-vol',self.vg,out) + + cfg = creators.astra_dict('BP'+self.appendString) + cfg['ProjectionDataId'] = sid + cfg['ReconstructionDataId'] = vid + cfg['ProjectorId'] = self.proj_id + bp_id = algorithm.create(cfg) + algorithm.run(bp_id) + + algorithm.delete(bp_id) + self.data_mod.delete([vid,sid]) + return out + + + + class OpTomoTranspose(scipy.sparse.linalg.LinearOperator): """This object provides the transpose operation (``.T``) of the OpTomo object. -- cgit v1.2.3 From ed717202a0c917958892e26322d6ea5173f7b32c Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 25 Apr 2016 17:04:39 +0200 Subject: Use FP/BP out argument in sample plugin --- samples/python/s018_plugin.py | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/samples/python/s018_plugin.py b/samples/python/s018_plugin.py index 31cca95..85b5486 100644 --- a/samples/python/s018_plugin.py +++ b/samples/python/s018_plugin.py @@ -30,30 +30,38 @@ import six # Define the plugin class (has to subclass astra.plugin.base) # Note that usually, these will be defined in a separate package/module -class SIRTPlugin(astra.plugin.base): - """Example of an ASTRA plugin class, implementing a simple 2D SIRT algorithm. +class LandweberPlugin(astra.plugin.base): + """Example of an ASTRA plugin class, implementing a simple 2D Landweber algorithm. Options: - 'rel_factor': relaxation factor (optional) + 'Relaxation': relaxation factor (optional) """ # The astra_name variable defines the name to use to # call the plugin from ASTRA - astra_name = "SIRT-PLUGIN" + astra_name = "LANDWEBER-PLUGIN" - def initialize(self,cfg, rel_factor = 1): + def initialize(self,cfg, Relaxation = 1): self.W = astra.OpTomo(cfg['ProjectorId']) self.vid = cfg['ReconstructionDataId'] self.sid = cfg['ProjectionDataId'] - self.rel = rel_factor + self.rel = Relaxation def run(self, its): v = astra.data2d.get_shared(self.vid) s = astra.data2d.get_shared(self.sid) + tv = np.zeros(v.shape, dtype=np.float32) + ts = np.zeros(s.shape, dtype=np.float32) W = self.W for i in range(its): - v[:] += self.rel*(W.T*(s - (W*v).reshape(s.shape))).reshape(v.shape)/s.size + W.FP(v,out=ts) + ts -= s # ts = W*v - s + + W.BP(ts,out=tv) + tv *= self.rel / s.size + + v -= tv # v = v - rel * W'*(W*v-s) / s.size if __name__=='__main__': @@ -75,20 +83,20 @@ if __name__=='__main__': # First we import the package that contains the plugin import s018_plugin # Then, we register the plugin class with ASTRA - astra.plugin.register(s018_plugin.SIRTPlugin) + astra.plugin.register(s018_plugin.LandweberPlugin) # Get a list of registered plugins six.print_(astra.plugin.get_registered()) # To get help on a registered plugin, use get_help - six.print_(astra.plugin.get_help('SIRT-PLUGIN')) + six.print_(astra.plugin.get_help('LANDWEBER-PLUGIN')) # Create data structures sid = astra.data2d.create('-sino', proj_geom, sinogram) vid = astra.data2d.create('-vol', vol_geom) # Create config using plugin name - cfg = astra.astra_dict('SIRT-PLUGIN') + cfg = astra.astra_dict('LANDWEBER-PLUGIN') cfg['ProjectorId'] = proj_id cfg['ProjectionDataId'] = sid cfg['ReconstructionDataId'] = vid @@ -103,18 +111,18 @@ if __name__=='__main__': rec = astra.data2d.get(vid) # Options for the plugin go in cfg['option'] - cfg = astra.astra_dict('SIRT-PLUGIN') + cfg = astra.astra_dict('LANDWEBER-PLUGIN') cfg['ProjectorId'] = proj_id cfg['ProjectionDataId'] = sid cfg['ReconstructionDataId'] = vid cfg['option'] = {} - cfg['option']['rel_factor'] = 1.5 + cfg['option']['Relaxation'] = 1.5 alg_id_rel = astra.algorithm.create(cfg) astra.algorithm.run(alg_id_rel, 100) rec_rel = astra.data2d.get(vid) # We can also use OpTomo to call the plugin - rec_op = W.reconstruct('SIRT-PLUGIN', sinogram, 100, extraOptions={'rel_factor':1.5}) + rec_op = W.reconstruct('LANDWEBER-PLUGIN', sinogram, 100, extraOptions={'Relaxation':1.5}) import pylab as pl pl.gray() -- cgit v1.2.3 From 73ad6a97531b6bc1d311eceb6ba2770fdf407b81 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 17 May 2016 11:22:17 +0200 Subject: Avoid defining singleton construction in header file That way, the call to new is always executed by code inside libastra. This avoids the situation where a singleton gets created by a copy of the constructor linked into an object file outside of libastra, such as a .mex file, which would then also cause the vtable to be outside of libastra. This situation would cause issues when .mex files are unloaded. --- include/astra/Singleton.h | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/include/astra/Singleton.h b/include/astra/Singleton.h index a256187..1ef4bba 100644 --- a/include/astra/Singleton.h +++ b/include/astra/Singleton.h @@ -57,15 +57,17 @@ class Singleton { m_singleton = 0; } + static void construct(); + // get singleton static T& getSingleton() { if (!m_singleton) - m_singleton = new T(); + construct(); return *m_singleton; } static T* getSingletonPtr() { if (!m_singleton) - m_singleton = new T(); + construct(); return m_singleton; } @@ -76,11 +78,23 @@ class Singleton { }; -#define DEFINE_SINGLETON(T) template<> T* Singleton::m_singleton = 0 +// We specifically avoid defining construct() in the header. +// That way, the call to new is always executed by code inside libastra. +// This avoids the situation where a singleton gets created by a copy +// of the constructor linked into an object file outside of libastra, such +// as a .mex file, which would then also cause the vtable to be outside of +// libastra. This situation would cause issues when .mex files are unloaded. + +#define DEFINE_SINGLETON(T) \ +template<> void Singleton::construct() { m_singleton = new T(); } \ +template<> T* Singleton::m_singleton = 0 + // This is a hack to support statements like // DEFINE_SINGLETON2(CTemplatedClass); -#define DEFINE_SINGLETON2(A,B) template<> A,B* Singleton::m_singleton = 0 +#define DEFINE_SINGLETON2(A,B) \ +template<> void Singleton::construct() { m_singleton = new A,B(); } \ +template<> A,B* Singleton::m_singleton = 0 } // end namespace -- cgit v1.2.3 From f9b68bafd90941d9faf53e5e2771361e3ab4336a Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 17 May 2016 16:00:28 +0200 Subject: Add sanity check --- include/astra/Singleton.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/astra/Singleton.h b/include/astra/Singleton.h index 1ef4bba..784c521 100644 --- a/include/astra/Singleton.h +++ b/include/astra/Singleton.h @@ -86,14 +86,14 @@ class Singleton { // libastra. This situation would cause issues when .mex files are unloaded. #define DEFINE_SINGLETON(T) \ -template<> void Singleton::construct() { m_singleton = new T(); } \ +template<> void Singleton::construct() { assert(!m_singleton); m_singleton = new T(); } \ template<> T* Singleton::m_singleton = 0 // This is a hack to support statements like // DEFINE_SINGLETON2(CTemplatedClass); #define DEFINE_SINGLETON2(A,B) \ -template<> void Singleton::construct() { m_singleton = new A,B(); } \ +template<> void Singleton::construct() { assert(!m_singleton); m_singleton = new A,B(); } \ template<> A,B* Singleton::m_singleton = 0 } // end namespace -- cgit v1.2.3 From b85245ff01ba2595b09f6bc606bb5676c129a499 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 20 Jun 2016 12:08:33 +0200 Subject: Improve volume block reduction The previous version would make the blocks too large due to inefficient computation of overlap. --- src/CompositeGeometryManager.cpp | 297 ++++++++++++++++++++++----------------- 1 file changed, 165 insertions(+), 132 deletions(-) diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index 084ba8c..eba6762 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -196,6 +196,69 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div return true; } + +static std::pair reduceProjectionVertical(const CVolumeGeometry3D* pVolGeom, const CProjectionGeometry3D* pProjGeom) +{ + double vmin_g, vmax_g; + + // reduce self to only cover intersection with projection of VolumePart + // (Project corners of volume, take bounding box) + + for (int i = 0; i < pProjGeom->getProjectionCount(); ++i) { + + double vol_u[8]; + double vol_v[8]; + + double pixx = pVolGeom->getPixelLengthX(); + double pixy = pVolGeom->getPixelLengthY(); + double pixz = pVolGeom->getPixelLengthZ(); + + // TODO: Is 0.5 sufficient? + double xmin = pVolGeom->getWindowMinX() - 0.5 * pixx; + double xmax = pVolGeom->getWindowMaxX() + 0.5 * pixx; + double ymin = pVolGeom->getWindowMinY() - 0.5 * pixy; + double ymax = pVolGeom->getWindowMaxY() + 0.5 * pixy; + double zmin = pVolGeom->getWindowMinZ() - 0.5 * pixz; + double zmax = pVolGeom->getWindowMaxZ() + 0.5 * pixz; + + pProjGeom->projectPoint(xmin, ymin, zmin, i, vol_u[0], vol_v[0]); + pProjGeom->projectPoint(xmin, ymin, zmax, i, vol_u[1], vol_v[1]); + pProjGeom->projectPoint(xmin, ymax, zmin, i, vol_u[2], vol_v[2]); + pProjGeom->projectPoint(xmin, ymax, zmax, i, vol_u[3], vol_v[3]); + pProjGeom->projectPoint(xmax, ymin, zmin, i, vol_u[4], vol_v[4]); + pProjGeom->projectPoint(xmax, ymin, zmax, i, vol_u[5], vol_v[5]); + pProjGeom->projectPoint(xmax, ymax, zmin, i, vol_u[6], vol_v[6]); + pProjGeom->projectPoint(xmax, ymax, zmax, i, vol_u[7], vol_v[7]); + + double vmin = vol_v[0]; + double vmax = vol_v[0]; + + for (int j = 1; j < 8; ++j) { + if (vol_v[j] < vmin) + vmin = vol_v[j]; + if (vol_v[j] > vmax) + vmax = vol_v[j]; + } + + if (i == 0 || vmin < vmin_g) + vmin_g = vmin; + if (i == 0 || vmax > vmax_g) + vmax_g = vmax; + } + + if (vmin_g < -1.0) + vmin_g = -1.0; + if (vmax_g > pProjGeom->getDetectorRowCount()) + vmax_g = pProjGeom->getDetectorRowCount(); + + if (vmin_g >= vmax_g) + vmin_g = vmax_g = 0.0; + + return std::pair(vmin_g, vmax_g); + +} + + CCompositeGeometryManager::CPart::CPart(const CPart& other) { eType = other.eType; @@ -236,119 +299,135 @@ size_t CCompositeGeometryManager::CPart::getSize() } +static bool testVolumeRange(const std::pair& fullRange, + const CVolumeGeometry3D *pVolGeom, + const CProjectionGeometry3D *pProjGeom, + int zmin, int zmax) +{ + double pixz = pVolGeom->getPixelLengthZ(); + + CVolumeGeometry3D test(pVolGeom->getGridColCount(), + pVolGeom->getGridRowCount(), + zmax - zmin, + pVolGeom->getWindowMinX(), + pVolGeom->getWindowMinY(), + pVolGeom->getWindowMinZ() + zmin * pixz, + pVolGeom->getWindowMaxX(), + pVolGeom->getWindowMaxY(), + pVolGeom->getWindowMinZ() + zmax * pixz); + + + std::pair subRange = reduceProjectionVertical(&test, pProjGeom); + + // empty + if (subRange.first == subRange.second) + return true; + + // fully outside of fullRange + if (subRange.first >= fullRange.second || subRange.second <= fullRange.first) + return true; + + return false; +} + CCompositeGeometryManager::CPart* CCompositeGeometryManager::CVolumePart::reduce(const CPart *_other) { const CProjectionPart *other = dynamic_cast(_other); assert(other); - // TODO: Is 0.5 sufficient? - double umin = -0.5; - double umax = other->pGeom->getDetectorColCount() + 0.5; - double vmin = -0.5; - double vmax = other->pGeom->getDetectorRowCount() + 0.5; - - double uu[4]; - double vv[4]; - uu[0] = umin; vv[0] = vmin; - uu[1] = umin; vv[1] = vmax; - uu[2] = umax; vv[2] = vmin; - uu[3] = umax; vv[3] = vmax; - - double pixx = pGeom->getPixelLengthX(); - double pixy = pGeom->getPixelLengthY(); - double pixz = pGeom->getPixelLengthZ(); - double xmin = pGeom->getWindowMinX() - 0.5 * pixx; - double xmax = pGeom->getWindowMaxX() + 0.5 * pixx; - double ymin = pGeom->getWindowMinY() - 0.5 * pixy; - double ymax = pGeom->getWindowMaxY() + 0.5 * pixy; - - // NB: Flipped - double zmax = pGeom->getWindowMinZ() - 2.5 * pixz; - double zmin = pGeom->getWindowMaxZ() + 2.5 * pixz; - - // TODO: This isn't as tight as it could be. - // In particular it won't detect the detector being - // missed entirely on the u side. - - for (int i = 0; i < other->pGeom->getProjectionCount(); ++i) { - for (int j = 0; j < 4; ++j) { - double px, py, pz; - - other->pGeom->backprojectPointX(i, uu[j], vv[j], xmin, py, pz); - //ASTRA_DEBUG("%f %f (%f - %f)", py, pz, ymin, ymax); - if (pz < zmin) zmin = pz; - if (pz > zmax) zmax = pz; - other->pGeom->backprojectPointX(i, uu[j], vv[j], xmax, py, pz); - //ASTRA_DEBUG("%f %f (%f - %f)", py, pz, ymin, ymax); - if (pz < zmin) zmin = pz; - if (pz > zmax) zmax = pz; - - other->pGeom->backprojectPointY(i, uu[j], vv[j], ymin, px, pz); - //ASTRA_DEBUG("%f %f (%f - %f)", px, pz, xmin, xmax); - if (pz < zmin) zmin = pz; - if (pz > zmax) zmax = pz; - other->pGeom->backprojectPointY(i, uu[j], vv[j], ymax, px, pz); - //ASTRA_DEBUG("%f %f (%f - %f)", px, pz, xmin, xmax); - if (pz < zmin) zmin = pz; - if (pz > zmax) zmax = pz; + std::pair fullRange = reduceProjectionVertical(pGeom, other->pGeom); + + int top_slice = 0, bottom_slice = 0; + + if (fullRange.first < fullRange.second) { + + + // TOP SLICE + + int zmin = 0; + int zmax = pGeom->getGridSliceCount()-1; // (Don't try empty region) + + // Setting top slice to zmin is always valid. + + while (zmin < zmax) { + int zmid = (zmin + zmax + 1) / 2; + + bool ok = testVolumeRange(fullRange, pGeom, other->pGeom, + 0, zmid); + + ASTRA_DEBUG("binsearch min: [%d,%d], %d, %s", zmin, zmax, zmid, ok ? "ok" : "removed too much"); + + if (ok) + zmin = zmid; + else + zmax = zmid - 1; } - } - //ASTRA_DEBUG("coord extent: %f - %f", zmin, zmax); + top_slice = zmin; + - // Clip both zmin and zmax to get rid of extreme (or infinite) values - // NB: When individual pz values are +/-Inf, the sign is determined - // by ray direction and on which side of the face the ray passes. - if (zmin < pGeom->getWindowMinZ() - 2*pixz) - zmin = pGeom->getWindowMinZ() - 2*pixz; - if (zmin > pGeom->getWindowMaxZ() + 2*pixz) - zmin = pGeom->getWindowMaxZ() + 2*pixz; - if (zmax < pGeom->getWindowMinZ() - 2*pixz) - zmax = pGeom->getWindowMinZ() - 2*pixz; - if (zmax > pGeom->getWindowMaxZ() + 2*pixz) - zmax = pGeom->getWindowMaxZ() + 2*pixz; + // BOTTOM SLICE - zmin = (zmin - pixz - pGeom->getWindowMinZ()) / pixz; - zmax = (zmax + pixz - pGeom->getWindowMinZ()) / pixz; + zmin = top_slice + 1; // (Don't try empty region) + zmax = pGeom->getGridSliceCount(); - int _zmin = (int)floor(zmin); - int _zmax = (int)ceil(zmax); + // Setting bottom slice to zmax is always valid - //ASTRA_DEBUG("index extent: %d - %d", _zmin, _zmax); + while (zmin < zmax) { + int zmid = (zmin + zmax) / 2; - if (_zmin < 0) - _zmin = 0; - if (_zmax > pGeom->getGridSliceCount()) - _zmax = pGeom->getGridSliceCount(); + bool ok = testVolumeRange(fullRange, pGeom, other->pGeom, + zmid, pGeom->getGridSliceCount()); + + ASTRA_DEBUG("binsearch max: [%d,%d], %d, %s", zmin, zmax, zmid, ok ? "ok" : "removed too much"); + + if (ok) + zmax = zmid; + else + zmin = zmid + 1; + + } + + bottom_slice = zmax; - if (_zmax <= _zmin) { - _zmin = _zmax = 0; } - //ASTRA_DEBUG("adjusted extent: %d - %d", _zmin, _zmax); + + ASTRA_DEBUG("found extent: %d - %d", top_slice, bottom_slice); + + top_slice -= 1; + if (top_slice < 0) + top_slice = 0; + bottom_slice += 1; + if (bottom_slice >= pGeom->getGridSliceCount()) + bottom_slice = pGeom->getGridSliceCount(); + + ASTRA_DEBUG("adjusted extent: %d - %d", top_slice, bottom_slice); + + double pixz = pGeom->getPixelLengthZ(); CVolumePart *sub = new CVolumePart(); sub->subX = this->subX; sub->subY = this->subY; - sub->subZ = this->subZ + _zmin; + sub->subZ = this->subZ + top_slice; sub->pData = pData; - if (_zmin == _zmax) { + if (top_slice == bottom_slice) { sub->pGeom = 0; } else { sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(), pGeom->getGridRowCount(), - _zmax - _zmin, + bottom_slice - top_slice, pGeom->getWindowMinX(), pGeom->getWindowMinY(), - pGeom->getWindowMinZ() + _zmin * pixz, + pGeom->getWindowMinZ() + top_slice * pixz, pGeom->getWindowMaxX(), pGeom->getWindowMaxY(), - pGeom->getWindowMinZ() + _zmax * pixz); + pGeom->getWindowMinZ() + bottom_slice * pixz); } - ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d", this->subZ, this->subZ + pGeom->getGridSliceCount(), this->subZ + _zmin, this->subZ + _zmax); + ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d ( %f - %f )", this->subZ, this->subZ + pGeom->getGridSliceCount(), this->subZ + top_slice, this->subZ + bottom_slice, pGeom->getWindowMinZ() + top_slice * pixz, pGeom->getWindowMinZ() + bottom_slice * pixz); return sub; } @@ -742,62 +821,16 @@ void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, s } + CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::reduce(const CPart *_other) { const CVolumePart *other = dynamic_cast(_other); assert(other); - double vmin_g, vmax_g; - - // reduce self to only cover intersection with projection of VolumePart - // (Project corners of volume, take bounding box) - - for (int i = 0; i < pGeom->getProjectionCount(); ++i) { - - double vol_u[8]; - double vol_v[8]; - - double pixx = other->pGeom->getPixelLengthX(); - double pixy = other->pGeom->getPixelLengthY(); - double pixz = other->pGeom->getPixelLengthZ(); - - // TODO: Is 0.5 sufficient? - double xmin = other->pGeom->getWindowMinX() - 0.5 * pixx; - double xmax = other->pGeom->getWindowMaxX() + 0.5 * pixx; - double ymin = other->pGeom->getWindowMinY() - 0.5 * pixy; - double ymax = other->pGeom->getWindowMaxY() + 0.5 * pixy; - double zmin = other->pGeom->getWindowMinZ() - 0.5 * pixz; - double zmax = other->pGeom->getWindowMaxZ() + 0.5 * pixz; - - pGeom->projectPoint(xmin, ymin, zmin, i, vol_u[0], vol_v[0]); - pGeom->projectPoint(xmin, ymin, zmax, i, vol_u[1], vol_v[1]); - pGeom->projectPoint(xmin, ymax, zmin, i, vol_u[2], vol_v[2]); - pGeom->projectPoint(xmin, ymax, zmax, i, vol_u[3], vol_v[3]); - pGeom->projectPoint(xmax, ymin, zmin, i, vol_u[4], vol_v[4]); - pGeom->projectPoint(xmax, ymin, zmax, i, vol_u[5], vol_v[5]); - pGeom->projectPoint(xmax, ymax, zmin, i, vol_u[6], vol_v[6]); - pGeom->projectPoint(xmax, ymax, zmax, i, vol_u[7], vol_v[7]); - - double vmin = vol_v[0]; - double vmax = vol_v[0]; - - for (int j = 1; j < 8; ++j) { - if (vol_v[j] < vmin) - vmin = vol_v[j]; - if (vol_v[j] > vmax) - vmax = vol_v[j]; - } - - if (i == 0 || vmin < vmin_g) - vmin_g = vmin; - if (i == 0 || vmax > vmax_g) - vmax_g = vmax; - } - - // fprintf(stderr, "v extent: %f %f\n", vmin_g, vmax_g); - - int _vmin = (int)floor(vmin_g - 1.0f); - int _vmax = (int)ceil(vmax_g + 1.0f); + std::pair r = reduceProjectionVertical(other->pGeom, pGeom); + // fprintf(stderr, "v extent: %f %f\n", r.first, r.second); + int _vmin = (int)floor(r.first - 1.0); + int _vmax = (int)ceil(r.second + 1.0); if (_vmin < 0) _vmin = 0; if (_vmax > pGeom->getDetectorRowCount()) -- cgit v1.2.3 From 085ae1be5e1cd456cfed9027d60f28f7cc9caf6e Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 20 Jun 2016 15:34:15 +0200 Subject: Fix inefficient block split logic --- src/CompositeGeometryManager.cpp | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index eba6762..c63faaa 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -662,7 +662,9 @@ void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::T size_t m = std::min(maxSize / sliceSize, maxDim); size_t blockSize = computeLinearSplit(m, div, sliceCount); - int rem = sliceCount % blockSize; + int rem = blockSize - (sliceCount % blockSize); + if (rem == blockSize) + rem = 0; ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); @@ -709,7 +711,9 @@ void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::T size_t m = std::min(maxSize / sliceSize, maxDim); size_t blockSize = computeLinearSplit(m, div, sliceCount); - int rem = sliceCount % blockSize; + int rem = blockSize - (sliceCount % blockSize); + if (rem == blockSize) + rem = 0; ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); @@ -756,7 +760,9 @@ void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::T size_t m = std::min(maxSize / sliceSize, maxDim); size_t blockSize = computeLinearSplit(m, div, sliceCount); - int rem = sliceCount % blockSize; + int rem = blockSize - (sliceCount % blockSize); + if (rem == blockSize) + rem = 0; ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); @@ -870,7 +876,11 @@ void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManage size_t m = std::min(maxSize / sliceSize, maxDim); size_t blockSize = computeLinearSplit(m, div, sliceCount); - int rem = sliceCount % blockSize; + int rem = blockSize - (sliceCount % blockSize); + if (rem == blockSize) + rem = 0; + + ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); for (int x = -(rem / 2); x < sliceCount; x += blockSize) { int newsubX = x; @@ -912,7 +922,11 @@ void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManage size_t m = std::min(maxSize / sliceSize, maxDim); size_t blockSize = computeLinearSplit(m, div, sliceCount); - int rem = sliceCount % blockSize; + int rem = blockSize - (sliceCount % blockSize); + if (rem == blockSize) + rem = 0; + + ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); for (int z = -(rem / 2); z < sliceCount; z += blockSize) { int newsubZ = z; -- cgit v1.2.3 From 558820f1421e79b32e542c133c83683e58df6d5e Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 24 Jun 2016 15:28:47 +0200 Subject: Compute FBP filter in spatial domain --- cuda/2d/fft.cu | 46 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 39 insertions(+), 7 deletions(-) diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu index 2bfd493..2d259a9 100644 --- a/cuda/2d/fft.cu +++ b/cuda/2d/fft.cu @@ -35,6 +35,7 @@ $Id$ #include #include "../../include/astra/Logging.h" +#include "astra/Fourier.h" using namespace astra; @@ -303,16 +304,48 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, float * pfFilt = new float[_iFFTFourierDetectorCount]; float * pfW = new float[_iFFTFourierDetectorCount]; + // We cache one Fourier transform for repeated FBP's of the same size + static float *pfData = 0; + static int iFilterCacheSize = 0; + + if (!pfData || iFilterCacheSize != _iFFTRealDetectorCount) { + // Compute filter in spatial domain + + delete[] pfData; + pfData = new float[2*_iFFTRealDetectorCount]; + int *ip = new int[int(2+sqrt(_iFFTRealDetectorCount)+1)]; + ip[0] = 0; + float32 *w = new float32[_iFFTRealDetectorCount/2]; + + for (int i = 0; i < _iFFTRealDetectorCount; ++i) { + pfData[2*i+1] = 0.0f; + + if (i & 1) { + int j = i; + if (2*j > _iFFTRealDetectorCount) + j = _iFFTRealDetectorCount - j; + float f = M_PI * j; + pfData[2*i] = -1 / (f*f); + } else { + pfData[2*i] = 0.0f; + } + } + + pfData[0] = 0.25f; + + cdft(2*_iFFTRealDetectorCount, -1, pfData, ip, w); + delete[] ip; + delete[] w; + + iFilterCacheSize = _iFFTRealDetectorCount; + } + for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount; - // filt = 2*( 0:(order/2) )./order; - pfFilt[iDetectorIndex] = 2.0f * fRelIndex; - //pfFilt[iDetectorIndex] = 1.0f; - - // w = 2*pi*(0:size(filt,2)-1)/order - pfW[iDetectorIndex] = 3.1415f * 2.0f * fRelIndex; + pfFilt[iDetectorIndex] = 2.0f * pfData[2*iDetectorIndex]; + pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex; } switch(_eFilter) @@ -866,5 +899,4 @@ void downloadDebugFilterReal(float * _pfHostSinogram, int _iProjectionCount, free(pfHostFilter); } - #endif -- cgit v1.2.3 From 302b0c7eacf46c7315f2390e270ad20f5025c342 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 28 Jul 2016 15:32:38 +0200 Subject: Convert one more string parsing function to stringToDouble --- matlab/mex/mexHelpFunctions.cpp | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/matlab/mex/mexHelpFunctions.cpp b/matlab/mex/mexHelpFunctions.cpp index 13c4ade..b13dde3 100644 --- a/matlab/mex/mexHelpFunctions.cpp +++ b/matlab/mex/mexHelpFunctions.cpp @@ -402,16 +402,13 @@ mxArray* stringToMxArray(std::string input) return pVector; } - // number - char* end; - double content = ::strtod(input.c_str(), &end); - bool isnumber = !*end; - if (isnumber) { - return mxCreateDoubleScalar(content); + try { + // number + return mxCreateDoubleScalar(StringUtil::stringToDouble(input)); + } catch (const StringUtil::bad_cast &) { + // string + return mxCreateString(input.c_str()); } - - // string - return mxCreateString(input.c_str()); } //----------------------------------------------------------------------------------------- // turn a c++ map into a matlab struct -- cgit v1.2.3 From 584fb584816aefca42518c9a6075ac2df814dac6 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 28 Jul 2016 15:32:50 +0200 Subject: Replace use of boost::split by own function --- include/astra/Utilities.h | 3 +++ matlab/mex/mexHelpFunctions.cpp | 13 ++++-------- python/astra/src/PythonPluginAlgorithm.cpp | 10 ++++----- src/Utilities.cpp | 34 ++++++++++++++++++------------ 4 files changed, 31 insertions(+), 29 deletions(-) diff --git a/include/astra/Utilities.h b/include/astra/Utilities.h index 8d7c44d..22adfe2 100644 --- a/include/astra/Utilities.h +++ b/include/astra/Utilities.h @@ -85,6 +85,9 @@ _AstraExport std::string doubleToString(double f); template _AstraExport std::string toString(T f); + +_AstraExport void splitString(std::vector &items, const std::string& s, const char *delim); + } diff --git a/matlab/mex/mexHelpFunctions.cpp b/matlab/mex/mexHelpFunctions.cpp index b13dde3..d957aea 100644 --- a/matlab/mex/mexHelpFunctions.cpp +++ b/matlab/mex/mexHelpFunctions.cpp @@ -33,11 +33,6 @@ $Id$ #include "mexHelpFunctions.h" #include "astra/Utilities.h" -#include -#include -#include -#include - using namespace std; using namespace astra; @@ -362,8 +357,8 @@ mxArray* stringToMxArray(std::string input) // split rows std::vector row_strings; std::vector col_strings; - boost::split(row_strings, input, boost::is_any_of(";")); - boost::split(col_strings, row_strings[0], boost::is_any_of(",")); + StringUtil::splitString(row_strings, input, ";"); + StringUtil::splitString(col_strings, row_strings[0], ","); // get dimensions int rows = row_strings.size(); @@ -375,7 +370,7 @@ mxArray* stringToMxArray(std::string input) // loop elements for (unsigned int row = 0; row < rows; row++) { - boost::split(col_strings, row_strings[row], boost::is_any_of(",")); + StringUtil::splitString(col_strings, row_strings[row], ","); // check size for (unsigned int col = 0; col < col_strings.size(); col++) { out[col*rows + row] = StringUtil::stringToFloat(col_strings[col]); @@ -389,7 +384,7 @@ mxArray* stringToMxArray(std::string input) // split std::vector items; - boost::split(items, input, boost::is_any_of(",")); + StringUtil::splitString(items, input, ","); // init matrix mxArray* pVector = mxCreateDoubleMatrix(1, items.size(), mxREAL); diff --git a/python/astra/src/PythonPluginAlgorithm.cpp b/python/astra/src/PythonPluginAlgorithm.cpp index 617c0f4..893db94 100644 --- a/python/astra/src/PythonPluginAlgorithm.cpp +++ b/python/astra/src/PythonPluginAlgorithm.cpp @@ -31,8 +31,6 @@ along with the ASTRA Toolbox. If not, see . #include "astra/Logging.h" #include "astra/Utilities.h" -#include -#include #include #include #include @@ -146,7 +144,7 @@ CPythonPluginAlgorithmFactory::~CPythonPluginAlgorithmFactory(){ PyObject * getClassFromString(std::string str){ std::vector items; - boost::split(items, str, boost::is_any_of(".")); + StringUtil::splitString(items, str, "."); PyObject *pyclass = PyImport_ImportModule(items[0].c_str()); if(pyclass==NULL){ logPythonError(); @@ -303,10 +301,10 @@ PyObject * pyStringFromString(std::string str){ PyObject* stringToPythonValue(std::string str){ if(str.find(";")!=std::string::npos){ std::vector rows, row; - boost::split(rows, str, boost::is_any_of(";")); + StringUtil::splitString(rows, str, ";"); PyObject *mat = PyList_New(rows.size()); for(unsigned int i=0; i vec; - boost::split(vec, str, boost::is_any_of(",")); + StringUtil::splitString(vec, str, ","); PyObject *veclist = PyList_New(vec.size()); for(unsigned int i=0;i -#include -#include - #include #include #include @@ -84,18 +80,16 @@ std::vector stringToDoubleVector(const std::string &s) template std::vector stringToVector(const std::string& s) { - // split - std::vector items; - boost::split(items, s, boost::is_any_of(",;")); - - // init list std::vector out; - out.resize(items.size()); + size_t current = 0; + size_t next; + do { + next = s.find_first_of(",;", current); + std::string t = s.substr(current, next - current); + out.push_back(stringTo(t)); + current = next + 1; + } while (next != std::string::npos); - // loop elements - for (unsigned int i = 0; i < items.size(); i++) { - out[i] = stringTo(items[i]); - } return out; } @@ -120,6 +114,18 @@ std::string doubleToString(double f) template<> std::string toString(float f) { return floatToString(f); } template<> std::string toString(double f) { return doubleToString(f); } +void splitString(std::vector &items, const std::string& s, + const char *delim) +{ + items.clear(); + size_t current = 0; + size_t next; + do { + next = s.find_first_of(",;", current); + items.push_back(s.substr(current, next - current)); + current = next + 1; + } while (next != std::string::npos); +} } -- cgit v1.2.3 From 53249b3ad63f0d08b9862a75602acf263d230d77 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 28 Jul 2016 16:33:37 +0200 Subject: Remove leftover code from Singleton --- include/astra/Singleton.h | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/include/astra/Singleton.h b/include/astra/Singleton.h index 784c521..9d3c088 100644 --- a/include/astra/Singleton.h +++ b/include/astra/Singleton.h @@ -45,11 +45,7 @@ class Singleton { public: // constructor - Singleton() { - assert(!m_singleton); - int offset = (uintptr_t)(T*)1 - (uintptr_t)(Singleton*)(T*)1; - m_singleton = (T*)((uintptr_t)this + offset); - }; + Singleton() { } // destructor virtual ~Singleton() { -- cgit v1.2.3