From d58f6b51493c7931cbc02feb3ffeb0eca6ea3a4e Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 5 Jul 2018 16:13:37 +0200 Subject: Refactor a few filter-related functions out of cuda code --- astra_vc14.vcxproj | 3 +- astra_vc14.vcxproj.filters | 9 +- build/linux/Makefile.in | 1 + build/msvc/gen.py | 3 +- cuda/2d/fbp.cu | 2 +- cuda/2d/fft.cu | 383 +--------------- cuda/3d/fdk.cu | 2 +- .../astra/CudaFilteredBackProjectionAlgorithm.h | 3 +- include/astra/Filters.h | 71 +++ include/astra/cuda/2d/astra.h | 1 - include/astra/cuda/2d/fbp.h | 2 +- include/astra/cuda/2d/fbp_filters.h | 61 --- include/astra/cuda/2d/fft.h | 8 +- src/CudaFilteredBackProjectionAlgorithm.cpp | 118 +---- src/Filters.cpp | 484 +++++++++++++++++++++ 15 files changed, 582 insertions(+), 569 deletions(-) create mode 100644 include/astra/Filters.h delete mode 100644 include/astra/cuda/2d/fbp_filters.h create mode 100644 src/Filters.cpp diff --git a/astra_vc14.vcxproj b/astra_vc14.vcxproj index 64c08ee..893305f 100644 --- a/astra_vc14.vcxproj +++ b/astra_vc14.vcxproj @@ -509,6 +509,7 @@ + @@ -596,6 +597,7 @@ + @@ -656,7 +658,6 @@ - diff --git a/astra_vc14.vcxproj.filters b/astra_vc14.vcxproj.filters index 3e9ff9a..fef6a8a 100644 --- a/astra_vc14.vcxproj.filters +++ b/astra_vc14.vcxproj.filters @@ -168,6 +168,9 @@ Global & Other\source + + Global & Other\source + Global & Other\source @@ -434,6 +437,9 @@ Global & Other\headers + + Global & Other\headers + Global & Other\headers @@ -638,9 +644,6 @@ CUDA\cuda headers - - CUDA\cuda headers - CUDA\cuda headers diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in index 1627a2e..0b90bd9 100644 --- a/build/linux/Makefile.in +++ b/build/linux/Makefile.in @@ -142,6 +142,7 @@ BASE_OBJECTS=\ src/FanFlatProjectionGeometry2D.lo \ src/FanFlatVecProjectionGeometry2D.lo \ src/FilteredBackProjectionAlgorithm.lo \ + src/Filters.lo \ src/Float32Data2D.lo \ src/Float32Data3D.lo \ src/Float32Data3DMemory.lo \ diff --git a/build/msvc/gen.py b/build/msvc/gen.py index fcc12d2..47e7440 100644 --- a/build/msvc/gen.py +++ b/build/msvc/gen.py @@ -214,6 +214,7 @@ P_astra["filters"]["Global & Other\\source"] = [ "src\\AstraObjectManager.cpp", "src\\CompositeGeometryManager.cpp", "src\\Config.cpp", +"src\\Filters.cpp", "src\\Fourier.cpp", "src\\Globals.cpp", "src\\Logging.cpp", @@ -292,7 +293,6 @@ P_astra["filters"]["CUDA\\cuda headers"] = [ "include\\astra\\cuda\\2d\\em.h", "include\\astra\\cuda\\2d\\fan_bp.h", "include\\astra\\cuda\\2d\\fan_fp.h", -"include\\astra\\cuda\\2d\\fbp_filters.h", "include\\astra\\cuda\\2d\\fbp.h", "include\\astra\\cuda\\2d\\fft.h", "include\\astra\\cuda\\2d\\par_bp.h", @@ -354,6 +354,7 @@ P_astra["filters"]["Global & Other\\headers"] = [ "include\\astra\\clog.h", "include\\astra\\CompositeGeometryManager.h", "include\\astra\\Config.h", +"include\\astra\\Filters.h", "include\\astra\\Fourier.h", "include\\astra\\Globals.h", "include\\astra\\Logging.h", diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu index 48fb7dc..7a8d2e9 100644 --- a/cuda/2d/fbp.cu +++ b/cuda/2d/fbp.cu @@ -130,7 +130,7 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* case astra::FILTER_FLATTOP: case astra::FILTER_PARZEN: { - genFilter(_eFilter, _fD, dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount, _fFilterParameter); + genCuFFTFilter(_eFilter, _fD, dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount, _fFilterParameter); uploadComplexArrayToDevice(dims.iProjAngles, iFreqBinCount, pHostFilter, (cufftComplex*)D_filter); break; diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu index bd8cab5..4454746 100644 --- a/cuda/2d/fft.cu +++ b/cuda/2d/fft.cu @@ -300,387 +300,13 @@ void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter, } } -void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, +void genCuFFTFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount, float _fParameter /* = -1.0f */) { - 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; - - pfFilt[iDetectorIndex] = 2.0f * pfData[2*iDetectorIndex]; - pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex; - } - - switch(_eFilter) - { - case FILTER_RAMLAK: - { - // do nothing - break; - } - case FILTER_SHEPPLOGAN: - { - // filt(2:end) = filt(2:end) .* (sin(w(2:end)/(2*d))./(w(2:end)/(2*d))) - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (sinf(pfW[iDetectorIndex] / 2.0f / _fD) / (pfW[iDetectorIndex] / 2.0f / _fD)); - } - break; - } - case FILTER_COSINE: - { - // filt(2:end) = filt(2:end) .* cos(w(2:end)/(2*d)) - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * cosf(pfW[iDetectorIndex] / 2.0f / _fD); - } - break; - } - case FILTER_HAMMING: - { - // filt(2:end) = filt(2:end) .* (.54 + .46 * cos(w(2:end)/d)) - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * ( 0.54f + 0.46f * cosf(pfW[iDetectorIndex] / _fD)); - } - break; - } - case FILTER_HANN: - { - // filt(2:end) = filt(2:end) .*(1+cos(w(2:end)./d)) / 2 - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (1.0f + cosf(pfW[iDetectorIndex] / _fD)) / 2.0f; - } - break; - } - case FILTER_TUKEY: - { - float fAlpha = _fParameter; - if(_fParameter < 0.0f) fAlpha = 0.5f; - float fN = (float)_iFFTFourierDetectorCount; - float fHalfN = fN / 2.0f; - float fEnumTerm = fAlpha * fHalfN; - float fDenom = (1.0f - fAlpha) * fHalfN; - float fBlockStart = fHalfN - fEnumTerm; - float fBlockEnd = fHalfN + fEnumTerm; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fAbsSmallN = fabs((float)iDetectorIndex); - float fStoredValue = 0.0f; - - if((fBlockStart <= fAbsSmallN) && (fAbsSmallN <= fBlockEnd)) - { - fStoredValue = 1.0f; - } - else - { - float fEnum = fAbsSmallN - fEnumTerm; - float fCosInput = M_PI * fEnum / fDenom; - fStoredValue = 0.5f * (1.0f + cosf(fCosInput)); - } - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_LANCZOS: - { - float fDenum = (float)(_iFFTFourierDetectorCount - 1); - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fX = 2.0f * fSmallN / fDenum - 1.0f; - float fSinInput = M_PI * fX; - float fStoredValue = 0.0f; - - if(fabsf(fSinInput) > 0.001f) - { - fStoredValue = sin(fSinInput)/fSinInput; - } - else - { - fStoredValue = 1.0f; - } - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_TRIANGULAR: - { - float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1); - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fAbsInput = fSmallN - fNMinusOne / 2.0f; - float fParenInput = fNMinusOne / 2.0f - fabsf(fAbsInput); - float fStoredValue = 2.0f / fNMinusOne * fParenInput; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_GAUSSIAN: - { - float fSigma = _fParameter; - if(_fParameter < 0.0f) fSigma = 0.4f; - float fN = (float)_iFFTFourierDetectorCount; - float fQuotient = (fN - 1.0f) / 2.0f; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fEnum = fSmallN - fQuotient; - float fDenom = fSigma * fQuotient; - float fPower = -0.5f * (fEnum / fDenom) * (fEnum / fDenom); - float fStoredValue = expf(fPower); - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_BARTLETTHANN: - { - const float fA0 = 0.62f; - const float fA1 = 0.48f; - const float fA2 = 0.38f; - float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fAbsInput = fSmallN / fNMinusOne - 0.5f; - float fFirstTerm = fA1 * fabsf(fAbsInput); - float fCosInput = 2.0f * M_PI * fSmallN / fNMinusOne; - float fSecondTerm = fA2 * cosf(fCosInput); - float fStoredValue = fA0 - fFirstTerm - fSecondTerm; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_BLACKMAN: - { - float fAlpha = _fParameter; - if(_fParameter < 0.0f) fAlpha = 0.16f; - float fA0 = (1.0f - fAlpha) / 2.0f; - float fA1 = 0.5f; - float fA2 = fAlpha / 2.0f; - float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1); - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fCosInput1 = 2.0f * M_PI * 0.5f * fSmallN / fNMinusOne; - float fCosInput2 = 4.0f * M_PI * 0.5f * fSmallN / fNMinusOne; - float fStoredValue = fA0 - fA1 * cosf(fCosInput1) + fA2 * cosf(fCosInput2); - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_NUTTALL: - { - const float fA0 = 0.355768f; - const float fA1 = 0.487396f; - const float fA2 = 0.144232f; - const float fA3 = 0.012604f; - float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fBaseCosInput = M_PI * fSmallN / fNMinusOne; - float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput); - float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput); - float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput); - float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_BLACKMANHARRIS: - { - const float fA0 = 0.35875f; - const float fA1 = 0.48829f; - const float fA2 = 0.14128f; - const float fA3 = 0.01168f; - float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fBaseCosInput = M_PI * fSmallN / fNMinusOne; - float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput); - float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput); - float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput); - float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_BLACKMANNUTTALL: - { - const float fA0 = 0.3635819f; - const float fA1 = 0.4891775f; - const float fA2 = 0.1365995f; - const float fA3 = 0.0106411f; - float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fBaseCosInput = M_PI * fSmallN / fNMinusOne; - float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput); - float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput); - float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput); - float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_FLATTOP: - { - const float fA0 = 1.0f; - const float fA1 = 1.93f; - const float fA2 = 1.29f; - const float fA3 = 0.388f; - const float fA4 = 0.032f; - float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fBaseCosInput = M_PI * fSmallN / fNMinusOne; - float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput); - float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput); - float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput); - float fFourthTerm = fA4 * cosf(8.0f * fBaseCosInput); - float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm + fFourthTerm; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_KAISER: - { - float fAlpha = _fParameter; - if(_fParameter < 0.0f) fAlpha = 3.0f; - float fPiTimesAlpha = M_PI * fAlpha; - float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1); - float fDenom = (float)j0((double)fPiTimesAlpha); - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fSquareInput = 2.0f * fSmallN / fNMinusOne - 1; - float fSqrtInput = 1.0f - fSquareInput * fSquareInput; - float fBesselInput = fPiTimesAlpha * sqrt(fSqrtInput); - float fEnum = (float)j0((double)fBesselInput); - float fStoredValue = fEnum / fDenom; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_PARZEN: - { - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fQ = fSmallN / (float)(_iFFTFourierDetectorCount - 1); - float fStoredValue = 0.0f; - - if(fQ <= 0.5f) - { - fStoredValue = 1.0f - 6.0f * fQ * fQ * (1.0f - fQ); - } - else - { - float fCubedValue = 1.0f - fQ; - fStoredValue = 2.0f * fCubedValue * fCubedValue * fCubedValue; - } - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - default: - { - ASTRA_ERROR("Cannot serve requested filter"); - } - } - - // filt(w>pi*d) = 0; - float fPiTimesD = M_PI * _fD; - for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fWValue = pfW[iDetectorIndex]; - - if(fWValue > fPiTimesD) - { - pfFilt[iDetectorIndex] = 0.0f; - } - } + float * pfFilt = astra::genFilter(_eFilter, _fD, _iProjectionCount, + _iFFTRealDetectorCount, + _iFFTFourierDetectorCount, _fParameter); for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { @@ -695,7 +321,6 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, } delete[] pfFilt; - delete[] pfW; } diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu index 46c07e7..194d2fb 100644 --- a/cuda/3d/fdk.cu +++ b/cuda/3d/fdk.cu @@ -253,7 +253,7 @@ bool FDK_Filter(cudaPitchedPtr D_projData, memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iHalfFFTSize); if (pfFilter == 0){ - astraCUDA::genFilter(astra::FILTER_RAMLAK, 1.0f, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize); + astraCUDA::genCuFFTFilter(astra::FILTER_RAMLAK, 1.0f, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize); } else { for (int i = 0; i < dims.iProjAngles * iHalfFFTSize; i++) { pHostFilter[i].x = pfFilter[i]; diff --git a/include/astra/CudaFilteredBackProjectionAlgorithm.h b/include/astra/CudaFilteredBackProjectionAlgorithm.h index 1280e9a..974914a 100644 --- a/include/astra/CudaFilteredBackProjectionAlgorithm.h +++ b/include/astra/CudaFilteredBackProjectionAlgorithm.h @@ -33,6 +33,7 @@ along with the ASTRA Toolbox. If not, see . #include "Float32ProjectionData2D.h" #include "Float32VolumeData2D.h" #include "CudaReconstructionAlgorithm2D.h" +#include "Filters.h" #include "cuda/2d/astra.h" @@ -52,8 +53,6 @@ private: float m_fFilterD; // frequency cut-off bool m_bShortScan; // short-scan mode for fan beam - static E_FBPFILTER _convertStringToFilter(const char * _filterType); - public: CCudaFilteredBackProjectionAlgorithm(); virtual ~CCudaFilteredBackProjectionAlgorithm(); diff --git a/include/astra/Filters.h b/include/astra/Filters.h new file mode 100644 index 0000000..ff3f5c8 --- /dev/null +++ b/include/astra/Filters.h @@ -0,0 +1,71 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2018, imec Vision Lab, University of Antwerp + 2014-2018, CWI, Amsterdam + +Contact: astra@astra-toolbox.com +Website: http://www.astra-toolbox.com/ + +This file is part of the ASTRA Toolbox. + + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +*/ + +#ifndef _INC_ASTRA_FILTERS_H +#define _INC_ASTRA_FILTERS_H + +namespace astra { + +enum E_FBPFILTER +{ + FILTER_ERROR, //< not a valid filter + FILTER_NONE, //< no filter (regular BP) + FILTER_RAMLAK, //< default FBP filter + FILTER_SHEPPLOGAN, //< Shepp-Logan + FILTER_COSINE, //< Cosine + FILTER_HAMMING, //< Hamming filter + FILTER_HANN, //< Hann filter + FILTER_TUKEY, //< Tukey filter + FILTER_LANCZOS, //< Lanczos filter + FILTER_TRIANGULAR, //< Triangular filter + FILTER_GAUSSIAN, //< Gaussian filter + FILTER_BARTLETTHANN, //< Bartlett-Hann filter + FILTER_BLACKMAN, //< Blackman filter + FILTER_NUTTALL, //< Nuttall filter, continuous first derivative + FILTER_BLACKMANHARRIS, //< Blackman-Harris filter + FILTER_BLACKMANNUTTALL, //< Blackman-Nuttall filter + FILTER_FLATTOP, //< Flat top filter + FILTER_KAISER, //< Kaiser filter + FILTER_PARZEN, //< Parzen filter + FILTER_PROJECTION, //< all projection directions share one filter + FILTER_SINOGRAM, //< every projection direction has its own filter + FILTER_RPROJECTION, //< projection filter in real space (as opposed to fourier space) + FILTER_RSINOGRAM, //< sinogram filter in real space + +}; + +// Generate filter of given size and parameters. Returns newly allocated array. +float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, + int _iFFTRealDetectorCount, + int _iFFTFourierDetectorCount, float _fParameter = -1.0f); + +// Convert string to filter type. Returns FILTER_ERROR if unrecognized. +E_FBPFILTER convertStringToFilter(const char * _filterType); + +} + +#endif diff --git a/include/astra/cuda/2d/astra.h b/include/astra/cuda/2d/astra.h index 6f0e2f0..a2f2219 100644 --- a/include/astra/cuda/2d/astra.h +++ b/include/astra/cuda/2d/astra.h @@ -28,7 +28,6 @@ along with the ASTRA Toolbox. If not, see . #ifndef _CUDA_ASTRA_H #define _CUDA_ASTRA_H -#include "fbp_filters.h" #include "dims.h" #include "algo.h" diff --git a/include/astra/cuda/2d/fbp.h b/include/astra/cuda/2d/fbp.h index 8666646..7febe69 100644 --- a/include/astra/cuda/2d/fbp.h +++ b/include/astra/cuda/2d/fbp.h @@ -26,7 +26,7 @@ along with the ASTRA Toolbox. If not, see . */ #include "algo.h" -#include "fbp_filters.h" +#include "astra/Filters.h" namespace astraCUDA { diff --git a/include/astra/cuda/2d/fbp_filters.h b/include/astra/cuda/2d/fbp_filters.h deleted file mode 100644 index 7c1121a..0000000 --- a/include/astra/cuda/2d/fbp_filters.h +++ /dev/null @@ -1,61 +0,0 @@ -/* ------------------------------------------------------------------------ -Copyright: 2010-2018, imec Vision Lab, University of Antwerp - 2014-2018, CWI, Amsterdam - -Contact: astra@astra-toolbox.com -Website: http://www.astra-toolbox.com/ - -This file is part of the ASTRA Toolbox. - - -The ASTRA Toolbox is free software: you can redistribute it and/or modify -it under the terms of the GNU General Public License as published by -the Free Software Foundation, either version 3 of the License, or -(at your option) any later version. - -The ASTRA Toolbox is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License -along with the ASTRA Toolbox. If not, see . - ------------------------------------------------------------------------ -*/ - -#ifndef FBP_FILTERS_H -#define FBP_FILTERS_H - -namespace astra { - -enum E_FBPFILTER -{ - FILTER_NONE, //< no filter (regular BP) - FILTER_RAMLAK, //< default FBP filter - FILTER_SHEPPLOGAN, //< Shepp-Logan - FILTER_COSINE, //< Cosine - FILTER_HAMMING, //< Hamming filter - FILTER_HANN, //< Hann filter - FILTER_TUKEY, //< Tukey filter - FILTER_LANCZOS, //< Lanczos filter - FILTER_TRIANGULAR, //< Triangular filter - FILTER_GAUSSIAN, //< Gaussian filter - FILTER_BARTLETTHANN, //< Bartlett-Hann filter - FILTER_BLACKMAN, //< Blackman filter - FILTER_NUTTALL, //< Nuttall filter, continuous first derivative - FILTER_BLACKMANHARRIS, //< Blackman-Harris filter - FILTER_BLACKMANNUTTALL, //< Blackman-Nuttall filter - FILTER_FLATTOP, //< Flat top filter - FILTER_KAISER, //< Kaiser filter - FILTER_PARZEN, //< Parzen filter - FILTER_PROJECTION, //< all projection directions share one filter - FILTER_SINOGRAM, //< every projection direction has its own filter - FILTER_RPROJECTION, //< projection filter in real space (as opposed to fourier space) - FILTER_RSINOGRAM, //< sinogram filter in real space -}; - -} - -#endif /* FBP_FILTERS_H */ diff --git a/include/astra/cuda/2d/fft.h b/include/astra/cuda/2d/fft.h index d36cae2..f5f8477 100644 --- a/include/astra/cuda/2d/fft.h +++ b/include/astra/cuda/2d/fft.h @@ -31,7 +31,7 @@ along with the ASTRA Toolbox. If not, see . #include #include -#include "fbp_filters.h" +#include "astra/Filters.h" namespace astraCUDA { @@ -60,9 +60,9 @@ void applyFilter(int _iProjectionCount, int _iFreqBinCount, int calcFFTFourierSize(int _iFFTRealSize); -void genFilter(astra::E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, - cufftComplex * _pFilter, int _iFFTRealDetectorCount, - int _iFFTFourierDetectorCount, float _fParameter = -1.0f); +void genCuFFTFilter(astra::E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, + cufftComplex * _pFilter, int _iFFTRealDetectorCount, + int _iFFTFourierDetectorCount, float _fParameter = -1.0f); void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount); diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp index 944f165..9ebbd60 100644 --- a/src/CudaFilteredBackProjectionAlgorithm.cpp +++ b/src/CudaFilteredBackProjectionAlgorithm.cpp @@ -27,10 +27,10 @@ along with the ASTRA Toolbox. If not, see . #include #include -#include #include "astra/AstraObjectManager.h" #include "astra/CudaProjector2D.h" +#include "astra/Filters.h" #include "astra/cuda/2d/astra.h" #include "astra/cuda/2d/fbp.h" @@ -75,7 +75,7 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) // filter type XMLNode node = _cfg.self.getSingleNode("FilterType"); if (node) - m_eFilter = _convertStringToFilter(node.getContent().c_str()); + m_eFilter = convertStringToFilter(node.getContent().c_str()); else m_eFilter = FILTER_RAMLAK; CC.markNodeParsed("FilterType"); @@ -215,6 +215,8 @@ bool CCudaFilteredBackProjectionAlgorithm::check() ASTRA_CONFIG_CHECK(m_pSinogram, "FBP_CUDA", "Invalid Projection Data Object."); ASTRA_CONFIG_CHECK(m_pReconstruction, "FBP_CUDA", "Invalid Reconstruction Data Object."); + ASTRA_CONFIG_CHECK(m_eFilter != FILTER_ERROR, "FBP_CUDA", "Invalid filter name."); + if((m_eFilter == FILTER_PROJECTION) || (m_eFilter == FILTER_SINOGRAM) || (m_eFilter == FILTER_RPROJECTION) || (m_eFilter == FILTER_RSINOGRAM)) { ASTRA_CONFIG_CHECK(m_pfFilter, "FBP_CUDA", "Invalid filter pointer."); @@ -235,116 +237,4 @@ bool CCudaFilteredBackProjectionAlgorithm::check() return true; } -static bool stringCompareLowerCase(const char * _stringA, const char * _stringB) -{ - int iCmpReturn = 0; - -#ifdef _MSC_VER - iCmpReturn = _stricmp(_stringA, _stringB); -#else - iCmpReturn = strcasecmp(_stringA, _stringB); -#endif - - return (iCmpReturn == 0); -} - -E_FBPFILTER CCudaFilteredBackProjectionAlgorithm::_convertStringToFilter(const char * _filterType) -{ - E_FBPFILTER output = FILTER_NONE; - - if(stringCompareLowerCase(_filterType, "ram-lak")) - { - output = FILTER_RAMLAK; - } - else if(stringCompareLowerCase(_filterType, "shepp-logan")) - { - output = FILTER_SHEPPLOGAN; - } - else if(stringCompareLowerCase(_filterType, "cosine")) - { - output = FILTER_COSINE; - } - else if(stringCompareLowerCase(_filterType, "hamming")) - { - output = FILTER_HAMMING; - } - else if(stringCompareLowerCase(_filterType, "hann")) - { - output = FILTER_HANN; - } - else if(stringCompareLowerCase(_filterType, "none")) - { - output = FILTER_NONE; - } - else if(stringCompareLowerCase(_filterType, "tukey")) - { - output = FILTER_TUKEY; - } - else if(stringCompareLowerCase(_filterType, "lanczos")) - { - output = FILTER_LANCZOS; - } - else if(stringCompareLowerCase(_filterType, "triangular")) - { - output = FILTER_TRIANGULAR; - } - else if(stringCompareLowerCase(_filterType, "gaussian")) - { - output = FILTER_GAUSSIAN; - } - else if(stringCompareLowerCase(_filterType, "barlett-hann")) - { - output = FILTER_BARTLETTHANN; - } - else if(stringCompareLowerCase(_filterType, "blackman")) - { - output = FILTER_BLACKMAN; - } - else if(stringCompareLowerCase(_filterType, "nuttall")) - { - output = FILTER_NUTTALL; - } - else if(stringCompareLowerCase(_filterType, "blackman-harris")) - { - output = FILTER_BLACKMANHARRIS; - } - else if(stringCompareLowerCase(_filterType, "blackman-nuttall")) - { - output = FILTER_BLACKMANNUTTALL; - } - else if(stringCompareLowerCase(_filterType, "flat-top")) - { - output = FILTER_FLATTOP; - } - else if(stringCompareLowerCase(_filterType, "kaiser")) - { - output = FILTER_KAISER; - } - else if(stringCompareLowerCase(_filterType, "parzen")) - { - output = FILTER_PARZEN; - } - else if(stringCompareLowerCase(_filterType, "projection")) - { - output = FILTER_PROJECTION; - } - else if(stringCompareLowerCase(_filterType, "sinogram")) - { - output = FILTER_SINOGRAM; - } - else if(stringCompareLowerCase(_filterType, "rprojection")) - { - output = FILTER_RPROJECTION; - } - else if(stringCompareLowerCase(_filterType, "rsinogram")) - { - output = FILTER_RSINOGRAM; - } - else - { - ASTRA_ERROR("Failed to convert \"%s\" into a filter.",_filterType); - } - - return output; -} diff --git a/src/Filters.cpp b/src/Filters.cpp new file mode 100644 index 0000000..30b2f24 --- /dev/null +++ b/src/Filters.cpp @@ -0,0 +1,484 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2018, imec Vision Lab, University of Antwerp + 2014-2018, CWI, Amsterdam + +Contact: astra@astra-toolbox.com +Website: http://www.astra-toolbox.com/ + +This file is part of the ASTRA Toolbox. + + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +*/ + +#include "astra/Globals.h" +#include "astra/Logging.h" +#include "astra/Fourier.h" +#include "astra/Filters.h" + +#include +#include + +namespace astra { + +float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, + int _iFFTRealDetectorCount, + int _iFFTFourierDetectorCount, float _fParameter /* = -1.0f */) +{ + 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; + + pfFilt[iDetectorIndex] = 2.0f * pfData[2*iDetectorIndex]; + pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex; + } + + switch(_eFilter) + { + case FILTER_RAMLAK: + { + // do nothing + break; + } + case FILTER_SHEPPLOGAN: + { + // filt(2:end) = filt(2:end) .* (sin(w(2:end)/(2*d))./(w(2:end)/(2*d))) + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (sinf(pfW[iDetectorIndex] / 2.0f / _fD) / (pfW[iDetectorIndex] / 2.0f / _fD)); + } + break; + } + case FILTER_COSINE: + { + // filt(2:end) = filt(2:end) .* cos(w(2:end)/(2*d)) + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * cosf(pfW[iDetectorIndex] / 2.0f / _fD); + } + break; + } + case FILTER_HAMMING: + { + // filt(2:end) = filt(2:end) .* (.54 + .46 * cos(w(2:end)/d)) + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * ( 0.54f + 0.46f * cosf(pfW[iDetectorIndex] / _fD)); + } + break; + } + case FILTER_HANN: + { + // filt(2:end) = filt(2:end) .*(1+cos(w(2:end)./d)) / 2 + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (1.0f + cosf(pfW[iDetectorIndex] / _fD)) / 2.0f; + } + break; + } + case FILTER_TUKEY: + { + float fAlpha = _fParameter; + if(_fParameter < 0.0f) fAlpha = 0.5f; + float fN = (float)_iFFTFourierDetectorCount; + float fHalfN = fN / 2.0f; + float fEnumTerm = fAlpha * fHalfN; + float fDenom = (1.0f - fAlpha) * fHalfN; + float fBlockStart = fHalfN - fEnumTerm; + float fBlockEnd = fHalfN + fEnumTerm; + + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fAbsSmallN = fabs((float)iDetectorIndex); + float fStoredValue = 0.0f; + + if((fBlockStart <= fAbsSmallN) && (fAbsSmallN <= fBlockEnd)) + { + fStoredValue = 1.0f; + } + else + { + float fEnum = fAbsSmallN - fEnumTerm; + float fCosInput = M_PI * fEnum / fDenom; + fStoredValue = 0.5f * (1.0f + cosf(fCosInput)); + } + + pfFilt[iDetectorIndex] *= fStoredValue; + } + + break; + } + case FILTER_LANCZOS: + { + float fDenum = (float)(_iFFTFourierDetectorCount - 1); + + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fSmallN = (float)iDetectorIndex; + float fX = 2.0f * fSmallN / fDenum - 1.0f; + float fSinInput = M_PI * fX; + float fStoredValue = 0.0f; + + if(fabsf(fSinInput) > 0.001f) + { + fStoredValue = sin(fSinInput)/fSinInput; + } + else + { + fStoredValue = 1.0f; + } + + pfFilt[iDetectorIndex] *= fStoredValue; + } + + break; + } + case FILTER_TRIANGULAR: + { + float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1); + + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fSmallN = (float)iDetectorIndex; + float fAbsInput = fSmallN - fNMinusOne / 2.0f; + float fParenInput = fNMinusOne / 2.0f - fabsf(fAbsInput); + float fStoredValue = 2.0f / fNMinusOne * fParenInput; + + pfFilt[iDetectorIndex] *= fStoredValue; + } + + break; + } + case FILTER_GAUSSIAN: + { + float fSigma = _fParameter; + if(_fParameter < 0.0f) fSigma = 0.4f; + float fN = (float)_iFFTFourierDetectorCount; + float fQuotient = (fN - 1.0f) / 2.0f; + + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fSmallN = (float)iDetectorIndex; + float fEnum = fSmallN - fQuotient; + float fDenom = fSigma * fQuotient; + float fPower = -0.5f * (fEnum / fDenom) * (fEnum / fDenom); + float fStoredValue = expf(fPower); + + pfFilt[iDetectorIndex] *= fStoredValue; + } + + break; + } + case FILTER_BARTLETTHANN: + { + const float fA0 = 0.62f; + const float fA1 = 0.48f; + const float fA2 = 0.38f; + float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; + + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fSmallN = (float)iDetectorIndex; + float fAbsInput = fSmallN / fNMinusOne - 0.5f; + float fFirstTerm = fA1 * fabsf(fAbsInput); + float fCosInput = 2.0f * M_PI * fSmallN / fNMinusOne; + float fSecondTerm = fA2 * cosf(fCosInput); + float fStoredValue = fA0 - fFirstTerm - fSecondTerm; + + pfFilt[iDetectorIndex] *= fStoredValue; + } + + break; + } + case FILTER_BLACKMAN: + { + float fAlpha = _fParameter; + if(_fParameter < 0.0f) fAlpha = 0.16f; + float fA0 = (1.0f - fAlpha) / 2.0f; + float fA1 = 0.5f; + float fA2 = fAlpha / 2.0f; + float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1); + + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fSmallN = (float)iDetectorIndex; + float fCosInput1 = 2.0f * M_PI * 0.5f * fSmallN / fNMinusOne; + float fCosInput2 = 4.0f * M_PI * 0.5f * fSmallN / fNMinusOne; + float fStoredValue = fA0 - fA1 * cosf(fCosInput1) + fA2 * cosf(fCosInput2); + + pfFilt[iDetectorIndex] *= fStoredValue; + } + + break; + } + case FILTER_NUTTALL: + { + const float fA0 = 0.355768f; + const float fA1 = 0.487396f; + const float fA2 = 0.144232f; + const float fA3 = 0.012604f; + float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; + + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fSmallN = (float)iDetectorIndex; + float fBaseCosInput = M_PI * fSmallN / fNMinusOne; + float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput); + float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput); + float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput); + float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm; + + pfFilt[iDetectorIndex] *= fStoredValue; + } + + break; + } + case FILTER_BLACKMANHARRIS: + { + const float fA0 = 0.35875f; + const float fA1 = 0.48829f; + const float fA2 = 0.14128f; + const float fA3 = 0.01168f; + float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; + + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fSmallN = (float)iDetectorIndex; + float fBaseCosInput = M_PI * fSmallN / fNMinusOne; + float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput); + float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput); + float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput); + float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm; + + pfFilt[iDetectorIndex] *= fStoredValue; + } + + break; + } + case FILTER_BLACKMANNUTTALL: + { + const float fA0 = 0.3635819f; + const float fA1 = 0.4891775f; + const float fA2 = 0.1365995f; + const float fA3 = 0.0106411f; + float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; + + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fSmallN = (float)iDetectorIndex; + float fBaseCosInput = M_PI * fSmallN / fNMinusOne; + float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput); + float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput); + float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput); + float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm; + + pfFilt[iDetectorIndex] *= fStoredValue; + } + + break; + } + case FILTER_FLATTOP: + { + const float fA0 = 1.0f; + const float fA1 = 1.93f; + const float fA2 = 1.29f; + const float fA3 = 0.388f; + const float fA4 = 0.032f; + float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; + + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fSmallN = (float)iDetectorIndex; + float fBaseCosInput = M_PI * fSmallN / fNMinusOne; + float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput); + float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput); + float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput); + float fFourthTerm = fA4 * cosf(8.0f * fBaseCosInput); + float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm + fFourthTerm; + + pfFilt[iDetectorIndex] *= fStoredValue; + } + + break; + } + case FILTER_KAISER: + { + float fAlpha = _fParameter; + if(_fParameter < 0.0f) fAlpha = 3.0f; + float fPiTimesAlpha = M_PI * fAlpha; + float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1); + float fDenom = (float)j0((double)fPiTimesAlpha); + + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fSmallN = (float)iDetectorIndex; + float fSquareInput = 2.0f * fSmallN / fNMinusOne - 1; + float fSqrtInput = 1.0f - fSquareInput * fSquareInput; + float fBesselInput = fPiTimesAlpha * sqrt(fSqrtInput); + float fEnum = (float)j0((double)fBesselInput); + float fStoredValue = fEnum / fDenom; + + pfFilt[iDetectorIndex] *= fStoredValue; + } + + break; + } + case FILTER_PARZEN: + { + for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fSmallN = (float)iDetectorIndex; + float fQ = fSmallN / (float)(_iFFTFourierDetectorCount - 1); + float fStoredValue = 0.0f; + + if(fQ <= 0.5f) + { + fStoredValue = 1.0f - 6.0f * fQ * fQ * (1.0f - fQ); + } + else + { + float fCubedValue = 1.0f - fQ; + fStoredValue = 2.0f * fCubedValue * fCubedValue * fCubedValue; + } + + pfFilt[iDetectorIndex] *= fStoredValue; + } + + break; + } + default: + { + ASTRA_ERROR("Cannot serve requested filter"); + } + } + + // filt(w>pi*d) = 0; + float fPiTimesD = M_PI * _fD; + for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fWValue = pfW[iDetectorIndex]; + + if(fWValue > fPiTimesD) + { + pfFilt[iDetectorIndex] = 0.0f; + } + } + + delete[] pfW; + + return pfFilt; +} + +static bool stringCompareLowerCase(const char * _stringA, const char * _stringB) +{ + int iCmpReturn = 0; + +#ifdef _MSC_VER + iCmpReturn = _stricmp(_stringA, _stringB); +#else + iCmpReturn = strcasecmp(_stringA, _stringB); +#endif + + return (iCmpReturn == 0); +} + +struct FilterNameMapEntry { + const char *m_name; + E_FBPFILTER m_type; +}; + +E_FBPFILTER convertStringToFilter(const char * _filterType) +{ + + static const FilterNameMapEntry map[] = { + { "ram-lak", FILTER_RAMLAK }, + { "shepp-logan", FILTER_SHEPPLOGAN }, + { "cosine", FILTER_COSINE }, + { "hamming", FILTER_HAMMING}, + { "hann", FILTER_HANN}, + { "tukey", FILTER_TUKEY }, + { "lanczos", FILTER_LANCZOS}, + { "triangular", FILTER_TRIANGULAR}, + { "gaussian", FILTER_GAUSSIAN}, + { "barlett-hann", FILTER_BARTLETTHANN }, + { "blackman", FILTER_BLACKMAN}, + { "nuttall", FILTER_NUTTALL }, + { "blackman-harris", FILTER_BLACKMANHARRIS }, + { "blackman-nuttall", FILTER_BLACKMANNUTTALL }, + { "flat-top", FILTER_FLATTOP }, + { "kaiser", FILTER_KAISER }, + { "parzen", FILTER_PARZEN }, + { "projection", FILTER_PROJECTION }, + { "sinogram", FILTER_SINOGRAM }, + { "rprojection", FILTER_RPROJECTION }, + { "rsinogram", FILTER_RSINOGRAM }, + { "none", FILTER_NONE }, + { 0, FILTER_ERROR } }; + + const FilterNameMapEntry *i; + + for (i = &map[0]; i->m_name; ++i) + if (stringCompareLowerCase(_filterType, i->m_name)) + return i->m_type; + + ASTRA_ERROR("Failed to convert \"%s\" into a filter.",_filterType); + + return FILTER_ERROR; +} + + + +} -- cgit v1.2.3 From d32df3bcd1910f56195e828a0f7fba8fc04b90ab Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 12 Jul 2018 12:11:52 +0200 Subject: Refactor filter config --- cuda/2d/fbp.cu | 32 ++++---- cuda/2d/fft.cu | 8 +- cuda/3d/fdk.cu | 4 +- .../astra/CudaFilteredBackProjectionAlgorithm.h | 6 +- include/astra/Filters.h | 21 ++++- include/astra/cuda/2d/fbp.h | 4 +- include/astra/cuda/2d/fft.h | 4 +- src/CudaFilteredBackProjectionAlgorithm.cpp | 83 ++++--------------- src/Filters.cpp | 96 ++++++++++++++++++---- 9 files changed, 140 insertions(+), 118 deletions(-) diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu index 7a8d2e9..1574ccc 100644 --- a/cuda/2d/fbp.cu +++ b/cuda/2d/fbp.cu @@ -88,7 +88,7 @@ bool FBP::init() return true; } -bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* = NULL */, int _iFilterWidth /* = 0 */, float _fD /* = 1.0f */, float _fFilterParameter /* = -1.0f */) +bool FBP::setFilter(const astra::SFilterConfig &_cfg) { if (D_filter) { @@ -96,7 +96,7 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* D_filter = 0; } - if (_eFilter == astra::FILTER_NONE) + if (_cfg.m_eType == astra::FILTER_NONE) return true; // leave D_filter set to 0 @@ -108,7 +108,7 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* allocateComplexOnDevice(dims.iProjAngles, iFreqBinCount, (cufftComplex**)&D_filter); - switch(_eFilter) + switch(_cfg.m_eType) { case astra::FILTER_NONE: // handled above @@ -130,7 +130,7 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* case astra::FILTER_FLATTOP: case astra::FILTER_PARZEN: { - genCuFFTFilter(_eFilter, _fD, dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount, _fFilterParameter); + genCuFFTFilter(_cfg, dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount); uploadComplexArrayToDevice(dims.iProjAngles, iFreqBinCount, pHostFilter, (cufftComplex*)D_filter); break; @@ -138,11 +138,11 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* case astra::FILTER_PROJECTION: { // make sure the offered filter has the correct size - assert(_iFilterWidth == iFreqBinCount); + assert(_cfg.m_iCustomFilterWidth == iFreqBinCount); for(int iFreqBinIndex = 0; iFreqBinIndex < iFreqBinCount; iFreqBinIndex++) { - float fValue = _pfHostFilter[iFreqBinIndex]; + float fValue = _cfg.m_pfCustomFilter[iFreqBinIndex]; for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++) { @@ -156,13 +156,13 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* case astra::FILTER_SINOGRAM: { // make sure the offered filter has the correct size - assert(_iFilterWidth == iFreqBinCount); + assert(_cfg.m_iCustomFilterWidth == iFreqBinCount); for(int iFreqBinIndex = 0; iFreqBinIndex < iFreqBinCount; iFreqBinIndex++) { for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++) { - float fValue = _pfHostFilter[iFreqBinIndex + iProjectionIndex * _iFilterWidth]; + float fValue = _cfg.m_pfCustomFilter[iFreqBinIndex + iProjectionIndex * _cfg.m_iCustomFilterWidth]; pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].x = fValue; pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].y = 0.0f; @@ -178,16 +178,16 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* float * pfHostRealFilter = new float[iRealFilterElementCount]; memset(pfHostRealFilter, 0, sizeof(float) * iRealFilterElementCount); - int iUsedFilterWidth = min(_iFilterWidth, iFFTRealDetCount); - int iStartFilterIndex = (_iFilterWidth - iUsedFilterWidth) / 2; + int iUsedFilterWidth = min(_cfg.m_iCustomFilterWidth, iFFTRealDetCount); + int iStartFilterIndex = (_cfg.m_iCustomFilterWidth - iUsedFilterWidth) / 2; int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth; - int iFilterShiftSize = _iFilterWidth / 2; + int iFilterShiftSize = _cfg.m_iCustomFilterWidth / 2; for(int iDetectorIndex = iStartFilterIndex; iDetectorIndex < iMaxFilterIndex; iDetectorIndex++) { int iFFTInFilterIndex = (iDetectorIndex + iFFTRealDetCount - iFilterShiftSize) % iFFTRealDetCount; - float fValue = _pfHostFilter[iDetectorIndex]; + float fValue = _cfg.m_pfCustomFilter[iDetectorIndex]; for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++) { @@ -213,11 +213,11 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* float* pfHostRealFilter = new float[iRealFilterElementCount]; memset(pfHostRealFilter, 0, sizeof(float) * iRealFilterElementCount); - int iUsedFilterWidth = min(_iFilterWidth, iFFTRealDetCount); - int iStartFilterIndex = (_iFilterWidth - iUsedFilterWidth) / 2; + int iUsedFilterWidth = min(_cfg.m_iCustomFilterWidth, iFFTRealDetCount); + int iStartFilterIndex = (_cfg.m_iCustomFilterWidth - iUsedFilterWidth) / 2; int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth; - int iFilterShiftSize = _iFilterWidth / 2; + int iFilterShiftSize = _cfg.m_iCustomFilterWidth / 2; for(int iDetectorIndex = iStartFilterIndex; iDetectorIndex < iMaxFilterIndex; iDetectorIndex++) { @@ -225,7 +225,7 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++) { - float fValue = _pfHostFilter[iDetectorIndex + iProjectionIndex * _iFilterWidth]; + float fValue = _cfg.m_pfCustomFilter[iDetectorIndex + iProjectionIndex * _cfg.m_iCustomFilterWidth]; pfHostRealFilter[iFFTInFilterIndex + iProjectionIndex * iFFTRealDetCount] = fValue; } } diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu index 4454746..864e325 100644 --- a/cuda/2d/fft.cu +++ b/cuda/2d/fft.cu @@ -300,13 +300,13 @@ void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter, } } -void genCuFFTFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, +void genCuFFTFilter(const SFilterConfig &_cfg, int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, - int _iFFTFourierDetectorCount, float _fParameter /* = -1.0f */) + int _iFFTFourierDetectorCount) { - float * pfFilt = astra::genFilter(_eFilter, _fD, _iProjectionCount, + float * pfFilt = astra::genFilter(_cfg, _iProjectionCount, _iFFTRealDetectorCount, - _iFFTFourierDetectorCount, _fParameter); + _iFFTFourierDetectorCount); for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu index 194d2fb..014529b 100644 --- a/cuda/3d/fdk.cu +++ b/cuda/3d/fdk.cu @@ -253,7 +253,9 @@ bool FDK_Filter(cudaPitchedPtr D_projData, memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iHalfFFTSize); if (pfFilter == 0){ - astraCUDA::genCuFFTFilter(astra::FILTER_RAMLAK, 1.0f, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize); + astra::SFilterConfig filter; + filter.m_eType = astra::FILTER_RAMLAK; + astraCUDA::genCuFFTFilter(filter, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize); } else { for (int i = 0; i < dims.iProjAngles * iHalfFFTSize; i++) { pHostFilter[i].x = pfFilter[i]; diff --git a/include/astra/CudaFilteredBackProjectionAlgorithm.h b/include/astra/CudaFilteredBackProjectionAlgorithm.h index 974914a..8ef5a57 100644 --- a/include/astra/CudaFilteredBackProjectionAlgorithm.h +++ b/include/astra/CudaFilteredBackProjectionAlgorithm.h @@ -46,11 +46,7 @@ public: static std::string type; private: - E_FBPFILTER m_eFilter; - float * m_pfFilter; - int m_iFilterWidth; // number of elements per projection direction in filter - float m_fFilterParameter; // some filters allow for parameterization (value < 0.0f -> no parameter) - float m_fFilterD; // frequency cut-off + SFilterConfig m_filterConfig; bool m_bShortScan; // short-scan mode for fan beam public: diff --git a/include/astra/Filters.h b/include/astra/Filters.h index ff3f5c8..eec2ba2 100644 --- a/include/astra/Filters.h +++ b/include/astra/Filters.h @@ -30,6 +30,9 @@ along with the ASTRA Toolbox. If not, see . namespace astra { +struct Config; +class CAlgorithm; + enum E_FBPFILTER { FILTER_ERROR, //< not a valid filter @@ -58,14 +61,28 @@ enum E_FBPFILTER }; +struct SFilterConfig { + E_FBPFILTER m_eType; + float m_fD; + float m_fParameter; + + float *m_pfCustomFilter; + int m_iCustomFilterWidth; + + SFilterConfig() : m_eType(FILTER_ERROR), m_fD(1.0f), m_fParameter(-1.0f), + m_pfCustomFilter(0), m_iCustomFilterWidth(0) { }; +}; + // Generate filter of given size and parameters. Returns newly allocated array. -float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, +float *genFilter(const SFilterConfig &_cfg, int _iProjectionCount, int _iFFTRealDetectorCount, - int _iFFTFourierDetectorCount, float _fParameter = -1.0f); + int _iFFTFourierDetectorCount); // Convert string to filter type. Returns FILTER_ERROR if unrecognized. E_FBPFILTER convertStringToFilter(const char * _filterType); +SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg); + } #endif diff --git a/include/astra/cuda/2d/fbp.h b/include/astra/cuda/2d/fbp.h index 7febe69..1adf3b1 100644 --- a/include/astra/cuda/2d/fbp.h +++ b/include/astra/cuda/2d/fbp.h @@ -75,9 +75,7 @@ public: // FILTER_COSINE, FILTER_HAMMING, and FILTER_HANN) // have a D variable, which gives the cutoff point in the frequency domain. // Setting this value to 1.0 will include the whole filter - bool setFilter(astra::E_FBPFILTER _eFilter, - const float * _pfHostFilter = NULL, - int _iFilterWidth = 0, float _fD = 1.0f, float _fFilterParameter = -1.0f); + bool setFilter(const astra::SFilterConfig &_cfg); bool setShortScan(bool ss) { m_bShortScan = ss; return true; } diff --git a/include/astra/cuda/2d/fft.h b/include/astra/cuda/2d/fft.h index f5f8477..33612a0 100644 --- a/include/astra/cuda/2d/fft.h +++ b/include/astra/cuda/2d/fft.h @@ -60,9 +60,9 @@ void applyFilter(int _iProjectionCount, int _iFreqBinCount, int calcFFTFourierSize(int _iFFTRealSize); -void genCuFFTFilter(astra::E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, +void genCuFFTFilter(const astra::SFilterConfig &_cfg, int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, - int _iFFTFourierDetectorCount, float _fParameter = -1.0f); + int _iFFTFourierDetectorCount); void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount); diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp index 9ebbd60..90b831e 100644 --- a/src/CudaFilteredBackProjectionAlgorithm.cpp +++ b/src/CudaFilteredBackProjectionAlgorithm.cpp @@ -45,15 +45,12 @@ CCudaFilteredBackProjectionAlgorithm::CCudaFilteredBackProjectionAlgorithm() { m_bIsInitialized = false; CCudaReconstructionAlgorithm2D::_clear(); - m_pfFilter = NULL; - m_fFilterParameter = -1.0f; - m_fFilterD = 1.0f; } CCudaFilteredBackProjectionAlgorithm::~CCudaFilteredBackProjectionAlgorithm() { - delete[] m_pfFilter; - m_pfFilter = NULL; + delete[] m_filterConfig.m_pfCustomFilter; + m_filterConfig.m_pfCustomFilter = NULL; } bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) @@ -71,59 +68,7 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) if (!m_bIsInitialized) return false; - - // filter type - XMLNode node = _cfg.self.getSingleNode("FilterType"); - if (node) - m_eFilter = convertStringToFilter(node.getContent().c_str()); - else - m_eFilter = FILTER_RAMLAK; - CC.markNodeParsed("FilterType"); - - // filter - node = _cfg.self.getSingleNode("FilterSinogramId"); - if (node) - { - int id = node.getContentInt(); - const CFloat32ProjectionData2D * pFilterData = dynamic_cast(CData2DManager::getSingleton().get(id)); - m_iFilterWidth = pFilterData->getGeometry()->getDetectorCount(); - int iFilterProjectionCount = pFilterData->getGeometry()->getProjectionAngleCount(); - - m_pfFilter = new float[m_iFilterWidth * iFilterProjectionCount]; - memcpy(m_pfFilter, pFilterData->getDataConst(), sizeof(float) * m_iFilterWidth * iFilterProjectionCount); - } - else - { - m_iFilterWidth = 0; - m_pfFilter = NULL; - } - CC.markNodeParsed("FilterSinogramId"); // TODO: Only for some types! - - // filter parameter - node = _cfg.self.getSingleNode("FilterParameter"); - if (node) - { - float fParameter = node.getContentNumerical(); - m_fFilterParameter = fParameter; - } - else - { - m_fFilterParameter = -1.0f; - } - CC.markNodeParsed("FilterParameter"); // TODO: Only for some types! - - // D value - node = _cfg.self.getSingleNode("FilterD"); - if (node) - { - float fD = node.getContentNumerical(); - m_fFilterD = fD; - } - else - { - m_fFilterD = 1.0f; - } - CC.markNodeParsed("FilterD"); // TODO: Only for some types! + m_filterConfig = getFilterConfigForAlgorithm(_cfg, this); // Fan beam short scan mode if (m_pSinogram && dynamic_cast(m_pSinogram->getGeometry())) { @@ -153,8 +98,8 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D * m_pReconstruction = _pReconstruction; m_iGPUIndex = _iGPUIndex; - m_eFilter = _eFilter; - m_iFilterWidth = _iFilterWidth; + m_filterConfig.m_eType = _eFilter; + m_filterConfig.m_iCustomFilterWidth = _iFilterWidth; m_bShortScan = false; // success @@ -167,7 +112,7 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D * { int iFilterElementCount = 0; - if((_eFilter != FILTER_SINOGRAM) && (_eFilter != FILTER_RSINOGRAM)) + if((m_filterConfig.m_eType != FILTER_SINOGRAM) && (m_filterConfig.m_eType != FILTER_RSINOGRAM)) { iFilterElementCount = _iFilterWidth; } @@ -176,15 +121,15 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D * iFilterElementCount = m_pSinogram->getAngleCount(); } - m_pfFilter = new float[iFilterElementCount]; - memcpy(m_pfFilter, _pfFilter, iFilterElementCount * sizeof(float)); + m_filterConfig.m_pfCustomFilter = new float[iFilterElementCount]; + memcpy(m_filterConfig.m_pfCustomFilter, _pfFilter, iFilterElementCount * sizeof(float)); } else { - m_pfFilter = NULL; + m_filterConfig.m_pfCustomFilter = NULL; } - m_fFilterParameter = _fFilterParameter; + m_filterConfig.m_fParameter = _fFilterParameter; return check(); } @@ -196,7 +141,7 @@ void CCudaFilteredBackProjectionAlgorithm::initCUDAAlgorithm() astraCUDA::FBP* pFBP = dynamic_cast(m_pAlgo); - bool ok = pFBP->setFilter(m_eFilter, m_pfFilter, m_iFilterWidth, m_fFilterD, m_fFilterParameter); + bool ok = pFBP->setFilter(m_filterConfig); if (!ok) { ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set filter"); ASTRA_ASSERT(ok); @@ -215,11 +160,11 @@ bool CCudaFilteredBackProjectionAlgorithm::check() ASTRA_CONFIG_CHECK(m_pSinogram, "FBP_CUDA", "Invalid Projection Data Object."); ASTRA_CONFIG_CHECK(m_pReconstruction, "FBP_CUDA", "Invalid Reconstruction Data Object."); - ASTRA_CONFIG_CHECK(m_eFilter != FILTER_ERROR, "FBP_CUDA", "Invalid filter name."); + ASTRA_CONFIG_CHECK(m_filterConfig.m_eType != FILTER_ERROR, "FBP_CUDA", "Invalid filter name."); - if((m_eFilter == FILTER_PROJECTION) || (m_eFilter == FILTER_SINOGRAM) || (m_eFilter == FILTER_RPROJECTION) || (m_eFilter == FILTER_RSINOGRAM)) + if((m_filterConfig.m_eType == FILTER_PROJECTION) || (m_filterConfig.m_eType == FILTER_SINOGRAM) || (m_filterConfig.m_eType == FILTER_RPROJECTION) || (m_filterConfig.m_eType == FILTER_RSINOGRAM)) { - ASTRA_CONFIG_CHECK(m_pfFilter, "FBP_CUDA", "Invalid filter pointer."); + ASTRA_CONFIG_CHECK(m_filterConfig.m_pfCustomFilter, "FBP_CUDA", "Invalid filter pointer."); } // check initializations diff --git a/src/Filters.cpp b/src/Filters.cpp index 30b2f24..bbd7cfd 100644 --- a/src/Filters.cpp +++ b/src/Filters.cpp @@ -29,15 +29,18 @@ along with the ASTRA Toolbox. If not, see . #include "astra/Logging.h" #include "astra/Fourier.h" #include "astra/Filters.h" +#include "astra/Config.h" +#include "astra/Float32ProjectionData2D.h" +#include "astra/AstraObjectManager.h" #include #include namespace astra { -float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, +float *genFilter(const SFilterConfig &_cfg, int _iProjectionCount, int _iFFTRealDetectorCount, - int _iFFTFourierDetectorCount, float _fParameter /* = -1.0f */) + int _iFFTFourierDetectorCount) { float * pfFilt = new float[_iFFTFourierDetectorCount]; float * pfW = new float[_iFFTFourierDetectorCount]; @@ -86,7 +89,7 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex; } - switch(_eFilter) + switch(_cfg.m_eType) { case FILTER_RAMLAK: { @@ -98,7 +101,7 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, // filt(2:end) = filt(2:end) .* (sin(w(2:end)/(2*d))./(w(2:end)/(2*d))) for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { - pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (sinf(pfW[iDetectorIndex] / 2.0f / _fD) / (pfW[iDetectorIndex] / 2.0f / _fD)); + pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (sinf(pfW[iDetectorIndex] / 2.0f / _cfg.m_fD) / (pfW[iDetectorIndex] / 2.0f / _cfg.m_fD)); } break; } @@ -107,7 +110,7 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, // filt(2:end) = filt(2:end) .* cos(w(2:end)/(2*d)) for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { - pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * cosf(pfW[iDetectorIndex] / 2.0f / _fD); + pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * cosf(pfW[iDetectorIndex] / 2.0f / _cfg.m_fD); } break; } @@ -116,7 +119,7 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, // filt(2:end) = filt(2:end) .* (.54 + .46 * cos(w(2:end)/d)) for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { - pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * ( 0.54f + 0.46f * cosf(pfW[iDetectorIndex] / _fD)); + pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * ( 0.54f + 0.46f * cosf(pfW[iDetectorIndex] / _cfg.m_fD)); } break; } @@ -125,14 +128,14 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, // filt(2:end) = filt(2:end) .*(1+cos(w(2:end)./d)) / 2 for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { - pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (1.0f + cosf(pfW[iDetectorIndex] / _fD)) / 2.0f; + pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (1.0f + cosf(pfW[iDetectorIndex] / _cfg.m_fD)) / 2.0f; } break; } case FILTER_TUKEY: { - float fAlpha = _fParameter; - if(_fParameter < 0.0f) fAlpha = 0.5f; + float fAlpha = _cfg.m_fParameter; + if(_cfg.m_fParameter < 0.0f) fAlpha = 0.5f; float fN = (float)_iFFTFourierDetectorCount; float fHalfN = fN / 2.0f; float fEnumTerm = fAlpha * fHalfN; @@ -204,8 +207,8 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, } case FILTER_GAUSSIAN: { - float fSigma = _fParameter; - if(_fParameter < 0.0f) fSigma = 0.4f; + float fSigma = _cfg.m_fParameter; + if(_cfg.m_fParameter < 0.0f) fSigma = 0.4f; float fN = (float)_iFFTFourierDetectorCount; float fQuotient = (fN - 1.0f) / 2.0f; @@ -245,8 +248,8 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, } case FILTER_BLACKMAN: { - float fAlpha = _fParameter; - if(_fParameter < 0.0f) fAlpha = 0.16f; + float fAlpha = _cfg.m_fParameter; + if(_cfg.m_fParameter < 0.0f) fAlpha = 0.16f; float fA0 = (1.0f - fAlpha) / 2.0f; float fA1 = 0.5f; float fA2 = fAlpha / 2.0f; @@ -356,8 +359,8 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, } case FILTER_KAISER: { - float fAlpha = _fParameter; - if(_fParameter < 0.0f) fAlpha = 3.0f; + float fAlpha = _cfg.m_fParameter; + if(_cfg.m_fParameter < 0.0f) fAlpha = 3.0f; float fPiTimesAlpha = M_PI * fAlpha; float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1); float fDenom = (float)j0((double)fPiTimesAlpha); @@ -406,7 +409,7 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, } // filt(w>pi*d) = 0; - float fPiTimesD = M_PI * _fD; + float fPiTimesD = M_PI * _cfg.m_fD; for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { float fWValue = pfW[iDetectorIndex]; @@ -480,5 +483,66 @@ E_FBPFILTER convertStringToFilter(const char * _filterType) } +SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg) +{ + ConfigStackCheck CC("getFilterConfig", _alg, _cfg); + + SFilterConfig c; + + // filter type + XMLNode node = _cfg.self.getSingleNode("FilterType"); + if (node) + c.m_eType = convertStringToFilter(node.getContent().c_str()); + else + c.m_eType = FILTER_RAMLAK; + CC.markNodeParsed("FilterType"); + + // filter + node = _cfg.self.getSingleNode("FilterSinogramId"); + if (node) + { + int id = node.getContentInt(); + const CFloat32ProjectionData2D * pFilterData = dynamic_cast(CData2DManager::getSingleton().get(id)); + c.m_iCustomFilterWidth = pFilterData->getGeometry()->getDetectorCount(); + int iFilterProjectionCount = pFilterData->getGeometry()->getProjectionAngleCount(); + + c.m_pfCustomFilter = new float[c.m_iCustomFilterWidth * iFilterProjectionCount]; + memcpy(c.m_pfCustomFilter, pFilterData->getDataConst(), sizeof(float) * c.m_iCustomFilterWidth * iFilterProjectionCount); + } + else + { + c.m_iCustomFilterWidth = 0; + c.m_pfCustomFilter = NULL; + } + CC.markNodeParsed("FilterSinogramId"); // TODO: Only for some types! + + // filter parameter + node = _cfg.self.getSingleNode("FilterParameter"); + if (node) + { + float fParameter = node.getContentNumerical(); + c.m_fParameter = fParameter; + } + else + { + c.m_fParameter = -1.0f; + } + CC.markNodeParsed("FilterParameter"); // TODO: Only for some types! + + // D value + node = _cfg.self.getSingleNode("FilterD"); + if (node) + { + float fD = node.getContentNumerical(); + c.m_fD = fD; + } + else + { + c.m_fD = 1.0f; + } + CC.markNodeParsed("FilterD"); // TODO: Only for some types! + + return c; +} } -- cgit v1.2.3 From da434892133fe0979c5e8ef8be277673452e7728 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 13 Jul 2018 14:48:19 +0200 Subject: Add filter size error reporting --- cuda/2d/fbp.cu | 7 +-- cuda/2d/fft.cu | 11 ----- cuda/3d/fdk.cu | 2 +- include/astra/Filters.h | 10 ++++- include/astra/cuda/2d/fft.h | 2 - src/CudaFDKAlgorithm3D.cpp | 4 +- src/CudaFilteredBackProjectionAlgorithm.cpp | 2 + src/Filters.cpp | 66 +++++++++++++++++++++++++++-- 8 files changed, 81 insertions(+), 23 deletions(-) diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu index 1574ccc..e8a224e 100644 --- a/cuda/2d/fbp.cu +++ b/cuda/2d/fbp.cu @@ -35,6 +35,7 @@ along with the ASTRA Toolbox. If not, see . #include "astra/cuda/3d/fdk.h" #include "astra/Logging.h" +#include "astra/Filters.h" #include @@ -55,7 +56,7 @@ static int calcNextPowerOfTwo(int n) int FBP::calcFourierFilterSize(int _iDetectorCount) { int iFFTRealDetCount = calcNextPowerOfTwo(2 * _iDetectorCount); - int iFreqBinCount = calcFFTFourierSize(iFFTRealDetCount); + int iFreqBinCount = astra::calcFFTFourierSize(iFFTRealDetCount); // CHECKME: Matlab makes this at least 64. Do we also need to? return iFreqBinCount; @@ -101,7 +102,7 @@ bool FBP::setFilter(const astra::SFilterConfig &_cfg) int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets); - int iFreqBinCount = calcFFTFourierSize(iFFTRealDetCount); + int iFreqBinCount = astra::calcFFTFourierSize(iFFTRealDetCount); cufftComplex * pHostFilter = new cufftComplex[dims.iProjAngles * iFreqBinCount]; memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iFreqBinCount); @@ -311,7 +312,7 @@ bool FBP::iterate(unsigned int iterations) if (D_filter) { int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets); - int iFFTFourDetCount = calcFFTFourierSize(iFFTRealDetCount); + int iFFTFourDetCount = astra::calcFFTFourierSize(iFFTRealDetCount); cufftComplex * pDevComplexSinogram = NULL; diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu index 864e325..dc0edc3 100644 --- a/cuda/2d/fft.cu +++ b/cuda/2d/fft.cu @@ -275,17 +275,6 @@ bool runCudaIFFT(int _iProjectionCount, const cufftComplex* _pDevSourceComplex, return true; } - -// Because the input is real, the Fourier transform is symmetric. -// CUFFT only outputs the first half (ignoring the redundant second half), -// and expects the same as input for the IFFT. -int calcFFTFourierSize(int _iFFTRealSize) -{ - int iFFTFourierSize = _iFFTRealSize / 2 + 1; - - return iFFTFourierSize; -} - void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount) { diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu index 014529b..1294721 100644 --- a/cuda/3d/fdk.cu +++ b/cuda/3d/fdk.cu @@ -246,7 +246,7 @@ bool FDK_Filter(cudaPitchedPtr D_projData, // Generate filter // TODO: Check errors int iPaddedDetCount = calcNextPowerOfTwo(2 * dims.iProjU); - int iHalfFFTSize = astraCUDA::calcFFTFourierSize(iPaddedDetCount); + int iHalfFFTSize = astra::calcFFTFourierSize(iPaddedDetCount); cufftComplex *pHostFilter = new cufftComplex[dims.iProjAngles * iHalfFFTSize]; diff --git a/include/astra/Filters.h b/include/astra/Filters.h index eec2ba2..bf0fb48 100644 --- a/include/astra/Filters.h +++ b/include/astra/Filters.h @@ -32,6 +32,7 @@ namespace astra { struct Config; class CAlgorithm; +class CProjectionGeometry2D; enum E_FBPFILTER { @@ -68,9 +69,11 @@ struct SFilterConfig { float *m_pfCustomFilter; int m_iCustomFilterWidth; + int m_iCustomFilterHeight; SFilterConfig() : m_eType(FILTER_ERROR), m_fD(1.0f), m_fParameter(-1.0f), - m_pfCustomFilter(0), m_iCustomFilterWidth(0) { }; + m_pfCustomFilter(0), m_iCustomFilterWidth(0), + m_iCustomFilterHeight(0) { }; }; // Generate filter of given size and parameters. Returns newly allocated array. @@ -83,6 +86,11 @@ E_FBPFILTER convertStringToFilter(const char * _filterType); SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg); +bool checkCustomFilterSize(const SFilterConfig &_cfg, const CProjectionGeometry2D &_geom); + +int calcFFTFourierSize(int _iFFTRealSize); + + } #endif diff --git a/include/astra/cuda/2d/fft.h b/include/astra/cuda/2d/fft.h index 33612a0..f77cbde 100644 --- a/include/astra/cuda/2d/fft.h +++ b/include/astra/cuda/2d/fft.h @@ -58,8 +58,6 @@ bool runCudaIFFT(int _iProjectionCount, const cufftComplex* _pDevSourceComplex, void applyFilter(int _iProjectionCount, int _iFreqBinCount, cufftComplex * _pSinogram, cufftComplex * _pFilter); -int calcFFTFourierSize(int _iFFTRealSize); - void genCuFFTFilter(const astra::SFilterConfig &_cfg, int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount); diff --git a/src/CudaFDKAlgorithm3D.cpp b/src/CudaFDKAlgorithm3D.cpp index 4d8fc5a..24ed04f 100644 --- a/src/CudaFDKAlgorithm3D.cpp +++ b/src/CudaFDKAlgorithm3D.cpp @@ -35,9 +35,9 @@ along with the ASTRA Toolbox. If not, see . #include "astra/CompositeGeometryManager.h" #include "astra/Logging.h" +#include "astra/Filters.h" #include "astra/cuda/3d/astra3d.h" -#include "astra/cuda/2d/fft.h" #include "astra/cuda/3d/util3d.h" using namespace std; @@ -156,7 +156,7 @@ bool CCudaFDKAlgorithm3D::initialize(const Config& _cfg) const CProjectionGeometry3D* projgeom = m_pSinogram->getGeometry(); const CProjectionGeometry2D* filtgeom = pFilterData->getGeometry(); int iPaddedDetCount = calcNextPowerOfTwo(2 * projgeom->getDetectorColCount()); - int iHalfFFTSize = astraCUDA::calcFFTFourierSize(iPaddedDetCount); + int iHalfFFTSize = calcFFTFourierSize(iPaddedDetCount); if(filtgeom->getDetectorCount()!=iHalfFFTSize || filtgeom->getProjectionAngleCount()!=projgeom->getProjectionCount()){ ASTRA_ERROR("Filter size does not match required size (%i angles, %i detectors)",projgeom->getProjectionCount(),iHalfFFTSize); return false; diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp index 90b831e..88e235b 100644 --- a/src/CudaFilteredBackProjectionAlgorithm.cpp +++ b/src/CudaFilteredBackProjectionAlgorithm.cpp @@ -176,6 +176,8 @@ bool CCudaFilteredBackProjectionAlgorithm::check() // check pixel supersampling ASTRA_CONFIG_CHECK(m_iPixelSuperSampling >= 0, "FBP_CUDA", "PixelSuperSampling must be a non-negative integer."); + ASTRA_CONFIG_CHECK(checkCustomFilterSize(m_filterConfig, *m_pSinogram->getGeometry()), "FBP_CUDA", "Filter size mismatch"); + // success m_bIsInitialized = true; diff --git a/src/Filters.cpp b/src/Filters.cpp index bbd7cfd..3ee3749 100644 --- a/src/Filters.cpp +++ b/src/Filters.cpp @@ -504,14 +504,15 @@ SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg) int id = node.getContentInt(); const CFloat32ProjectionData2D * pFilterData = dynamic_cast(CData2DManager::getSingleton().get(id)); c.m_iCustomFilterWidth = pFilterData->getGeometry()->getDetectorCount(); - int iFilterProjectionCount = pFilterData->getGeometry()->getProjectionAngleCount(); + c.m_iCustomFilterHeight = pFilterData->getGeometry()->getProjectionAngleCount(); - c.m_pfCustomFilter = new float[c.m_iCustomFilterWidth * iFilterProjectionCount]; - memcpy(c.m_pfCustomFilter, pFilterData->getDataConst(), sizeof(float) * c.m_iCustomFilterWidth * iFilterProjectionCount); + c.m_pfCustomFilter = new float[c.m_iCustomFilterWidth * c.m_iCustomFilterHeight]; + memcpy(c.m_pfCustomFilter, pFilterData->getDataConst(), sizeof(float) * c.m_iCustomFilterWidth * c.m_iCustomFilterHeight); } else { c.m_iCustomFilterWidth = 0; + c.m_iCustomFilterHeight = 0; c.m_pfCustomFilter = NULL; } CC.markNodeParsed("FilterSinogramId"); // TODO: Only for some types! @@ -545,4 +546,63 @@ SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg) return c; } +static int calcNextPowerOfTwo(int n) +{ + int x = 1; + while (x < n) + x *= 2; + + return x; +} + +// Because the input is real, the Fourier transform is symmetric. +// CUFFT only outputs the first half (ignoring the redundant second half), +// and expects the same as input for the IFFT. +int calcFFTFourierSize(int _iFFTRealSize) +{ + int iFFTFourierSize = _iFFTRealSize / 2 + 1; + + return iFFTFourierSize; +} + +bool checkCustomFilterSize(const SFilterConfig &_cfg, const CProjectionGeometry2D &_geom) { + int iExpectedWidth = -1, iExpectedHeight = 1; + + switch (_cfg.m_eType) { + case FILTER_ERROR: + ASTRA_ERROR("checkCustomFilterSize: internal error; FILTER_ERROR passed"); + return false; + case FILTER_NONE: + return true; + case FILTER_SINOGRAM: + iExpectedHeight = _geom.getProjectionAngleCount(); + // fallthrough + case FILTER_PROJECTION: + { + int iPaddedDetCount = calcNextPowerOfTwo(2 * _geom.getDetectorCount()); + iExpectedWidth = calcFFTFourierSize(iPaddedDetCount); + } + if (_cfg.m_iCustomFilterWidth != iExpectedWidth || + _cfg.m_iCustomFilterHeight != iExpectedHeight) + { + ASTRA_ERROR("filter size mismatch: %dx%d (received) is not %dx%d (expected)", _cfg.m_iCustomFilterHeight, _cfg.m_iCustomFilterWidth, iExpectedHeight, iExpectedWidth); + return false; + } + return true; + case FILTER_RSINOGRAM: + iExpectedHeight = _geom.getProjectionAngleCount(); + // fallthrough + case FILTER_RPROJECTION: + if (_cfg.m_iCustomFilterHeight != iExpectedHeight) + { + ASTRA_ERROR("filter size mismatch: %dx%d (received) is not %dxX (expected)", _cfg.m_iCustomFilterHeight, _cfg.m_iCustomFilterWidth, iExpectedHeight); + return false; + } + return true; + default: + // Non-custom filters; nothing to check. + return true; + } +} + } -- cgit v1.2.3 From bc65a0395e6c3930bac0440cc894990dd04cb704 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 13 Jul 2018 16:20:07 +0200 Subject: Reorganize more filter size functions --- cuda/2d/fbp.cu | 16 +++------------- cuda/2d/fft.cu | 2 +- include/astra/Filters.h | 4 +++- src/Filters.cpp | 6 +++--- 4 files changed, 10 insertions(+), 18 deletions(-) diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu index e8a224e..a5b8a7a 100644 --- a/cuda/2d/fbp.cu +++ b/cuda/2d/fbp.cu @@ -42,20 +42,10 @@ along with the ASTRA Toolbox. If not, see . namespace astraCUDA { - -static int calcNextPowerOfTwo(int n) -{ - int x = 1; - while (x < n) - x *= 2; - - return x; -} - // static int FBP::calcFourierFilterSize(int _iDetectorCount) { - int iFFTRealDetCount = calcNextPowerOfTwo(2 * _iDetectorCount); + int iFFTRealDetCount = astra::calcNextPowerOfTwo(2 * _iDetectorCount); int iFreqBinCount = astra::calcFFTFourierSize(iFFTRealDetCount); // CHECKME: Matlab makes this at least 64. Do we also need to? @@ -101,7 +91,7 @@ bool FBP::setFilter(const astra::SFilterConfig &_cfg) return true; // leave D_filter set to 0 - int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets); + int iFFTRealDetCount = astra::calcNextPowerOfTwo(2 * dims.iProjDets); int iFreqBinCount = astra::calcFFTFourierSize(iFFTRealDetCount); cufftComplex * pHostFilter = new cufftComplex[dims.iProjAngles * iFreqBinCount]; @@ -311,7 +301,7 @@ bool FBP::iterate(unsigned int iterations) if (D_filter) { - int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets); + int iFFTRealDetCount = astra::calcNextPowerOfTwo(2 * dims.iProjDets); int iFFTFourDetCount = astra::calcFFTFourierSize(iFFTRealDetCount); cufftComplex * pDevComplexSinogram = NULL; diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu index dc0edc3..2e94b79 100644 --- a/cuda/2d/fft.cu +++ b/cuda/2d/fft.cu @@ -293,7 +293,7 @@ void genCuFFTFilter(const SFilterConfig &_cfg, int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount) { - float * pfFilt = astra::genFilter(_cfg, _iProjectionCount, + float * pfFilt = astra::genFilter(_cfg, _iFFTRealDetectorCount, _iFFTFourierDetectorCount); diff --git a/include/astra/Filters.h b/include/astra/Filters.h index bf0fb48..a1dec97 100644 --- a/include/astra/Filters.h +++ b/include/astra/Filters.h @@ -77,7 +77,7 @@ struct SFilterConfig { }; // Generate filter of given size and parameters. Returns newly allocated array. -float *genFilter(const SFilterConfig &_cfg, int _iProjectionCount, +float *genFilter(const SFilterConfig &_cfg, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount); @@ -88,6 +88,8 @@ SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg); bool checkCustomFilterSize(const SFilterConfig &_cfg, const CProjectionGeometry2D &_geom); + +int calcNextPowerOfTwo(int _iValue); int calcFFTFourierSize(int _iFFTRealSize); diff --git a/src/Filters.cpp b/src/Filters.cpp index 3ee3749..c13aa6b 100644 --- a/src/Filters.cpp +++ b/src/Filters.cpp @@ -38,7 +38,7 @@ along with the ASTRA Toolbox. If not, see . namespace astra { -float *genFilter(const SFilterConfig &_cfg, int _iProjectionCount, +float *genFilter(const SFilterConfig &_cfg, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount) { @@ -546,10 +546,10 @@ SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg) return c; } -static int calcNextPowerOfTwo(int n) +int calcNextPowerOfTwo(int n) { int x = 1; - while (x < n) + while (x < n && x > 0) x *= 2; return x; -- cgit v1.2.3 From fa83c8a76fa19862de5c68914a5ef81397ae89a2 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 13 Jul 2018 16:21:01 +0200 Subject: Move CPU FBP to common filter code --- include/astra/FilteredBackProjectionAlgorithm.h | 5 ++ src/FilteredBackProjectionAlgorithm.cpp | 69 ++++++++++++++++++------- 2 files changed, 54 insertions(+), 20 deletions(-) diff --git a/include/astra/FilteredBackProjectionAlgorithm.h b/include/astra/FilteredBackProjectionAlgorithm.h index 1cd4296..a234845 100644 --- a/include/astra/FilteredBackProjectionAlgorithm.h +++ b/include/astra/FilteredBackProjectionAlgorithm.h @@ -35,6 +35,7 @@ along with the ASTRA Toolbox. If not, see . #include "Projector2D.h" #include "Float32ProjectionData2D.h" #include "Float32VolumeData2D.h" +#include "Filters.h" namespace astra { @@ -144,6 +145,10 @@ public: */ virtual std::string description() const; +protected: + + SFilterConfig m_filterConfig; + }; // inline functions diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp index bb2e722..ed72aa6 100644 --- a/src/FilteredBackProjectionAlgorithm.cpp +++ b/src/FilteredBackProjectionAlgorithm.cpp @@ -81,6 +81,9 @@ void CFilteredBackProjectionAlgorithm::clear() m_pSinogram = NULL; m_pReconstruction = NULL; m_bIsInitialized = false; + + delete[] m_filterConfig.m_pfCustomFilter; + m_filterConfig.m_pfCustomFilter = 0; } @@ -151,6 +154,9 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) delete[] projectionAngles; } + m_filterConfig = getFilterConfigForAlgorithm(_cfg, this); + + // TODO: check that the angles are linearly spaced between 0 and pi // success @@ -195,11 +201,14 @@ bool CFilteredBackProjectionAlgorithm::initialize(CProjector2D* _pProjector, CFloat32VolumeData2D* _pVolume, CFloat32ProjectionData2D* _pSinogram) { + clear(); + // store classes m_pProjector = _pProjector; m_pReconstruction = _pVolume; m_pSinogram = _pSinogram; + m_filterConfig = SFilterConfig(); // TODO: check that the angles are linearly spaced between 0 and pi @@ -214,6 +223,15 @@ bool CFilteredBackProjectionAlgorithm::_check() { ASTRA_CONFIG_CHECK(CReconstructionAlgorithm2D::_check(), "FBP", "Error in ReconstructionAlgorithm2D initialization"); + ASTRA_CONFIG_CHECK(m_filterConfig.m_eType != FILTER_ERROR, "FBP", "Invalid filter name."); + + if((m_filterConfig.m_eType == FILTER_PROJECTION) || (m_filterConfig.m_eType == FILTER_SINOGRAM) || (m_filterConfig.m_eType == FILTER_RPROJECTION) || (m_filterConfig.m_eType == FILTER_RSINOGRAM)) + { + ASTRA_CONFIG_CHECK(m_filterConfig.m_pfCustomFilter, "FBP", "Invalid filter pointer."); + } + + ASTRA_CONFIG_CHECK(checkCustomFilterSize(m_filterConfig, *m_pSinogram->getGeometry()), "FBP", "Filter size mismatch"); + // success return true; } @@ -249,29 +267,36 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D ASTRA_ASSERT(_pFilteredSinogram->getAngleCount() == m_pSinogram->getAngleCount()); ASTRA_ASSERT(_pFilteredSinogram->getDetectorCount() == m_pSinogram->getDetectorCount()); + ASTRA_ASSERT(m_filterConfig.m_eType != FILTER_ERROR); + if (m_filterConfig.m_eType == FILTER_NONE) + return; int iAngleCount = m_pProjector->getProjectionGeometry()->getProjectionAngleCount(); int iDetectorCount = m_pProjector->getProjectionGeometry()->getDetectorCount(); - // We'll zero-pad to the smallest power of two at least 64 and - // at least 2*iDetectorCount - int zpDetector = 64; - int nextPow2 = 5; - while (zpDetector < iDetectorCount*2) { - zpDetector *= 2; - nextPow2++; - } + int zpDetector = calcNextPowerOfTwo(2 * m_pSinogram->getDetectorCount()); + int iHalfFFTSize = astra::calcFFTFourierSize(zpDetector); // Create filter - float32* filter = new float32[zpDetector]; - - for (int iDetector = 0; iDetector <= zpDetector/2; iDetector++) - filter[iDetector] = (2.0f * iDetector)/zpDetector; - - for (int iDetector = zpDetector/2+1; iDetector < zpDetector; iDetector++) - filter[iDetector] = (2.0f * (zpDetector - iDetector)) / zpDetector; - + float *pfFilter = 0; + switch (m_filterConfig.m_eType) { + case FILTER_ERROR: + case FILTER_NONE: + // Should have been handled before + ASTRA_ASSERT(false); + return; + case FILTER_SINOGRAM: + case FILTER_PROJECTION: + case FILTER_RSINOGRAM: + case FILTER_RPROJECTION: + // TODO + ASTRA_ERROR("Unimplemented filter type"); + return; + + default: + pfFilter = genFilter(m_filterConfig, zpDetector, iHalfFFTSize); + } float32* pf = new float32[2 * iAngleCount * zpDetector]; int *ip = new int[int(2+sqrt((float)zpDetector)+1)]; @@ -301,9 +326,13 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D // Filter for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) { float32* pfRow = pf + iAngle * 2 * zpDetector; - for (int iDetector = 0; iDetector < zpDetector; ++iDetector) { - pfRow[2*iDetector] *= filter[iDetector]; - pfRow[2*iDetector+1] *= filter[iDetector]; + for (int iDetector = 0; iDetector < iHalfFFTSize; ++iDetector) { + pfRow[2*iDetector] *= pfFilter[iDetector]; + pfRow[2*iDetector+1] *= pfFilter[iDetector]; + } + for (int iDetector = iHalfFFTSize; iDetector < zpDetector; ++iDetector) { + pfRow[2*iDetector] *= pfFilter[zpDetector - iDetector]; + pfRow[2*iDetector+1] *= pfFilter[zpDetector - iDetector]; } } @@ -324,7 +353,7 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D delete[] pf; delete[] w; delete[] ip; - delete[] filter; + delete[] pfFilter; } } -- cgit v1.2.3 From db631e6cd00303c04ba5541ddf98a90d588a10c4 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 13 Jul 2018 16:31:12 +0200 Subject: Add SINOGRAM/PROJECTION filter modes to CPU FBP --- src/FilteredBackProjectionAlgorithm.cpp | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp index ed72aa6..49494d0 100644 --- a/src/FilteredBackProjectionAlgorithm.cpp +++ b/src/FilteredBackProjectionAlgorithm.cpp @@ -279,15 +279,22 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D int iHalfFFTSize = astra::calcFFTFourierSize(zpDetector); // Create filter + bool bFilterMultiAngle = false; float *pfFilter = 0; + float *pfFilter_delete = 0; switch (m_filterConfig.m_eType) { case FILTER_ERROR: case FILTER_NONE: // Should have been handled before ASTRA_ASSERT(false); return; - case FILTER_SINOGRAM: case FILTER_PROJECTION: + pfFilter = m_filterConfig.m_pfCustomFilter; + break; + case FILTER_SINOGRAM: + bFilterMultiAngle = true; + pfFilter = m_filterConfig.m_pfCustomFilter; + break; case FILTER_RSINOGRAM: case FILTER_RPROJECTION: // TODO @@ -296,6 +303,7 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D default: pfFilter = genFilter(m_filterConfig, zpDetector, iHalfFFTSize); + pfFilter_delete = pfFilter; } float32* pf = new float32[2 * iAngleCount * zpDetector]; @@ -326,13 +334,16 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D // Filter for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) { float32* pfRow = pf + iAngle * 2 * zpDetector; + float *pfFilterRow = pfFilter; + if (bFilterMultiAngle) + pfFilterRow += iAngle * iHalfFFTSize; for (int iDetector = 0; iDetector < iHalfFFTSize; ++iDetector) { - pfRow[2*iDetector] *= pfFilter[iDetector]; - pfRow[2*iDetector+1] *= pfFilter[iDetector]; + pfRow[2*iDetector] *= pfFilterRow[iDetector]; + pfRow[2*iDetector+1] *= pfFilterRow[iDetector]; } for (int iDetector = iHalfFFTSize; iDetector < zpDetector; ++iDetector) { - pfRow[2*iDetector] *= pfFilter[zpDetector - iDetector]; - pfRow[2*iDetector+1] *= pfFilter[zpDetector - iDetector]; + pfRow[2*iDetector] *= pfFilterRow[zpDetector - iDetector]; + pfRow[2*iDetector+1] *= pfFilterRow[zpDetector - iDetector]; } } @@ -353,7 +364,7 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D delete[] pf; delete[] w; delete[] ip; - delete[] pfFilter; + delete[] pfFilter_delete; } } -- cgit v1.2.3 From 9bce55d46758bc79ef2504f68cda6e79c81f4cba Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 13 Jul 2018 17:39:51 +0200 Subject: Add RSINOGRAM/RPROJECTION filter modes to CPU FBP --- src/FilteredBackProjectionAlgorithm.cpp | 85 ++++++++++++++++++++++++++------- 1 file changed, 68 insertions(+), 17 deletions(-) diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp index 49494d0..67a12a2 100644 --- a/src/FilteredBackProjectionAlgorithm.cpp +++ b/src/FilteredBackProjectionAlgorithm.cpp @@ -278,8 +278,14 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D int zpDetector = calcNextPowerOfTwo(2 * m_pSinogram->getDetectorCount()); int iHalfFFTSize = astra::calcFFTFourierSize(zpDetector); + // cdft setup + int *ip = new int[int(2+sqrt((float)zpDetector)+1)]; + ip[0] = 0; + float32 *w = new float32[zpDetector/2]; + // Create filter bool bFilterMultiAngle = false; + bool bFilterComplex = false; float *pfFilter = 0; float *pfFilter_delete = 0; switch (m_filterConfig.m_eType) { @@ -289,6 +295,8 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D ASTRA_ASSERT(false); return; case FILTER_PROJECTION: + // Fourier space, real, half the coefficients (because symmetric) + // 1 x iHalfFFTSize pfFilter = m_filterConfig.m_pfCustomFilter; break; case FILTER_SINOGRAM: @@ -296,20 +304,47 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D pfFilter = m_filterConfig.m_pfCustomFilter; break; case FILTER_RSINOGRAM: + bFilterMultiAngle = true; + // fall-through case FILTER_RPROJECTION: - // TODO - ASTRA_ERROR("Unimplemented filter type"); - return; + { + bFilterComplex = true; + + int count = bFilterMultiAngle ? iAngleCount : 1; + // Spatial, real, full convolution kernel + // Center in center (or right-of-center for even sized.) + // I.e., 0 1 0 and 0 0 1 0 both correspond to the identity + + pfFilter = new float[2 * zpDetector * count]; + pfFilter_delete = pfFilter; + + int iUsedFilterWidth = min(m_filterConfig.m_iCustomFilterWidth, zpDetector); + int iStartFilterIndex = (m_filterConfig.m_iCustomFilterWidth - iUsedFilterWidth) / 2; + int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth; + + int iFilterShiftSize = m_filterConfig.m_iCustomFilterWidth / 2; + + for (int i = 0; i < count; ++i) { + float *rOut = pfFilter + i * 2 * zpDetector; + float *rIn = m_filterConfig.m_pfCustomFilter + i * m_filterConfig.m_iCustomFilterWidth; + memset(rOut, 0, sizeof(float) * 2 * zpDetector); + + for(int j = iStartFilterIndex; j < iMaxFilterIndex; j++) { + int iFFTInFilterIndex = (j + zpDetector - iFilterShiftSize) % zpDetector; + rOut[2 * iFFTInFilterIndex] = rIn[j]; + } + + cdft(2*zpDetector, -1, rOut, ip, w); + } + break; + } default: pfFilter = genFilter(m_filterConfig, zpDetector, iHalfFFTSize); pfFilter_delete = pfFilter; } float32* pf = new float32[2 * iAngleCount * zpDetector]; - int *ip = new int[int(2+sqrt((float)zpDetector)+1)]; - ip[0]=0; - float32 *w = new float32[zpDetector/2]; // Copy and zero-pad data for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) { @@ -332,18 +367,34 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D } // Filter - for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) { - float32* pfRow = pf + iAngle * 2 * zpDetector; - float *pfFilterRow = pfFilter; - if (bFilterMultiAngle) - pfFilterRow += iAngle * iHalfFFTSize; - for (int iDetector = 0; iDetector < iHalfFFTSize; ++iDetector) { - pfRow[2*iDetector] *= pfFilterRow[iDetector]; - pfRow[2*iDetector+1] *= pfFilterRow[iDetector]; + if (bFilterComplex) { + for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) { + float32* pfRow = pf + iAngle * 2 * zpDetector; + float *pfFilterRow = pfFilter; + if (bFilterMultiAngle) + pfFilterRow += iAngle * 2 * zpDetector; + + for (int i = 0; i < zpDetector; ++i) { + float re = pfRow[2*i] * pfFilterRow[2*i] - pfRow[2*i+1] * pfFilterRow[2*i+1]; + float im = pfRow[2*i] * pfFilterRow[2*i+1] + pfRow[2*i+1] * pfFilterRow[2*i]; + pfRow[2*i] = re; + pfRow[2*i+1] = im; + } } - for (int iDetector = iHalfFFTSize; iDetector < zpDetector; ++iDetector) { - pfRow[2*iDetector] *= pfFilterRow[zpDetector - iDetector]; - pfRow[2*iDetector+1] *= pfFilterRow[zpDetector - iDetector]; + } else { + for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) { + float32* pfRow = pf + iAngle * 2 * zpDetector; + float *pfFilterRow = pfFilter; + if (bFilterMultiAngle) + pfFilterRow += iAngle * iHalfFFTSize; + for (int iDetector = 0; iDetector < iHalfFFTSize; ++iDetector) { + pfRow[2*iDetector] *= pfFilterRow[iDetector]; + pfRow[2*iDetector+1] *= pfFilterRow[iDetector]; + } + for (int iDetector = iHalfFFTSize; iDetector < zpDetector; ++iDetector) { + pfRow[2*iDetector] *= pfFilterRow[zpDetector - iDetector]; + pfRow[2*iDetector+1] *= pfFilterRow[zpDetector - iDetector]; + } } } -- cgit v1.2.3 From 4d741fc8e6c7930f7a8e27f54c55e0ad4949ed07 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 17 Jul 2018 16:54:13 +0200 Subject: Add sample scripts --- samples/matlab/s023_FBP_filters.m | 96 ++++++++++++++++++++++++++++++ samples/python/s023_FBP_filters.py | 116 +++++++++++++++++++++++++++++++++++++ 2 files changed, 212 insertions(+) create mode 100644 samples/matlab/s023_FBP_filters.m create mode 100644 samples/python/s023_FBP_filters.py diff --git a/samples/matlab/s023_FBP_filters.m b/samples/matlab/s023_FBP_filters.m new file mode 100644 index 0000000..4abec7e --- /dev/null +++ b/samples/matlab/s023_FBP_filters.m @@ -0,0 +1,96 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2018, imec Vision Lab, University of Antwerp +% 2014-2018, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@astra-toolbox.com +% Website: http://www.astra-toolbox.com/ +% ----------------------------------------------------------------------- + + +% This sample script illustrates three ways of passing filters to FBP. +% They work with both the FBP (CPU) and the FBP_CUDA (GPU) algorithms. + +N = 256; + +vol_geom = astra_create_vol_geom(N, N); +proj_geom = astra_create_proj_geom('parallel', 1.0, N, linspace2(0,pi,180)); + +proj_id = astra_create_projector('strip', proj_geom, vol_geom); + +P = phantom(256); + +[sinogram_id, sinogram] = astra_create_sino(P, proj_id); + +rec_id = astra_mex_data2d('create', '-vol', vol_geom); + +cfg = astra_struct('FBP'); +cfg.ReconstructionDataId = rec_id; +cfg.ProjectionDataId = sinogram_id; +cfg.ProjectorId = proj_id; + + +% 1. Use a standard Ram-Lak filter +cfg.FilterType = 'ram-lak'; + +alg_id = astra_mex_algorithm('create', cfg); +astra_mex_algorithm('run', alg_id); +rec_RL = astra_mex_data2d('get', rec_id); +astra_mex_algorithm('delete', alg_id); + + +% 2. Define a filter in Fourier space +% This is assumed to be symmetric, and ASTRA therefore expects only half + +% The full filter size should be the smallest power of two that is at least +% twice the number of detector pixels. +fullFilterSize = 2*N; +kernel = [linspace2(0, 1, floor(fullFilterSize / 2)) linspace2(1, 0, ceil(fullFilterSize / 2))]; +halfFilterSize = floor(fullFilterSize / 2) + 1; +filter = kernel(1:halfFilterSize); + +filter_geom = astra_create_proj_geom('parallel', 1.0, halfFilterSize, [0]); +filter_id = astra_mex_data2d('create', '-sino', filter_geom, filter); + +cfg.FilterType = 'projection'; +cfg.FilterSinogramId = filter_id; + +alg_id = astra_mex_algorithm('create', cfg); +astra_mex_algorithm('run', alg_id); +rec_filter = astra_mex_data2d('get', rec_id); +astra_mex_algorithm('delete', alg_id); + +% 3. Define a (spatial) convolution kernel directly +% For a kernel of odd size 2*k+1, the central component is at kernel(k+1) +% For a kernel of even size 2*k, the central component is at kernel(k+1) + +kernel = zeros(1, N); +for i = 0:floor(N/4)-1 + f = pi * (2*i + 1); + val = -2.0 / (f * f); + kernel(floor(N/2) + 1 + (2*i+1)) = val; + kernel(floor(N/2) + 1 - (2*i+1)) = val; +end +kernel(floor(N/2)+1) = 0.5; + +kernel_geom = astra_create_proj_geom('parallel', 1.0, N, [0]); +kernel_id = astra_mex_data2d('create', '-sino', kernel_geom, kernel); + +cfg.FilterType = 'rprojection'; +cfg.FilterSinogramId = kernel_id; + +alg_id = astra_mex_algorithm('create', cfg); +astra_mex_algorithm('run', alg_id); +rec_kernel = astra_mex_data2d('get', rec_id); +astra_mex_algorithm('delete', alg_id); + +figure(1); imshow(P, []); +figure(2); imshow(rec_RL, []); +figure(3); imshow(rec_filter, []); +figure(4); imshow(rec_kernel, []); + + +astra_mex_data2d('delete', rec_id); +astra_mex_data2d('delete', sinogram_id); +astra_mex_projector('delete', proj_id); diff --git a/samples/python/s023_FBP_filters.py b/samples/python/s023_FBP_filters.py new file mode 100644 index 0000000..11518ac --- /dev/null +++ b/samples/python/s023_FBP_filters.py @@ -0,0 +1,116 @@ +# ----------------------------------------------------------------------- +# Copyright: 2010-2018, imec Vision Lab, University of Antwerp +# 2013-2018, CWI, Amsterdam +# +# Contact: astra@astra-toolbox.com +# Website: http://www.astra-toolbox.com/ +# +# This file is part of the ASTRA Toolbox. +# +# +# The ASTRA Toolbox is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# The ASTRA Toolbox is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see . +# +# ----------------------------------------------------------------------- + +import astra +import numpy as np +import scipy.io + + +# This sample script illustrates three ways of passing filters to FBP. +# They work with both the FBP (CPU) and the FBP_CUDA (GPU) algorithms. + + +N = 256 + +vol_geom = astra.create_vol_geom(N, N) +proj_geom = astra.create_proj_geom('parallel', 1.0, N, np.linspace(0,np.pi,180,False)) + +P = scipy.io.loadmat('phantom.mat')['phantom256'] + +proj_id = astra.create_projector('strip',proj_geom,vol_geom) +sinogram_id, sinogram = astra.create_sino(P, proj_id) + +rec_id = astra.data2d.create('-vol', vol_geom) +cfg = astra.astra_dict('FBP') +cfg['ReconstructionDataId'] = rec_id +cfg['ProjectionDataId'] = sinogram_id +cfg['ProjectorId'] = proj_id + + + +# 1. Use a standard Ram-Lak filter +cfg['FilterType'] = 'ram-lak' + +alg_id = astra.algorithm.create(cfg) +astra.algorithm.run(alg_id) +rec_RL = astra.data2d.get(rec_id) +astra.algorithm.delete(alg_id) + +# 2. Define a filter in Fourier space +# This is assumed to be symmetric, and ASTRA therefore expects only half + +# The full filter size should be the smallest power of two that is at least +# twice the number of detector pixels. +fullFilterSize = 2*N +kernel = np.append( np.linspace(0, 1, fullFilterSize//2, endpoint=False), np.linspace(1, 0, fullFilterSize//2, endpoint=False) ) +halfFilterSize = fullFilterSize // 2 + 1 +filter = np.reshape(kernel[0:halfFilterSize], (1, halfFilterSize)) + +filter_geom = astra.create_proj_geom('parallel', 1.0, halfFilterSize, [0]); +filter_id = astra.data2d.create('-sino', filter_geom, filter); + +cfg['FilterType'] = 'projection' +cfg['FilterSinogramId'] = filter_id +alg_id = astra.algorithm.create(cfg) +astra.algorithm.run(alg_id) +rec_filter = astra.data2d.get(rec_id) +astra.algorithm.delete(alg_id) + + +# 3. Define a (spatial) convolution kernel directly +# For a kernel of odd size 2*k+1, the central component is at kernel[k] +# For a kernel of even size 2*k, the central component is at kernel[k] +kernel = np.zeros((1, N)) +for i in range(0,N//4): + f = np.pi * (2*i + 1) + val = -2.0 / (f * f) + kernel[0, N//2 + (2*i+1)] = val + kernel[0, N//2 - (2*i+1)] = val +kernel[0, N//2] = 0.5 +kernel_geom = astra.create_proj_geom('parallel', 1.0, N, [0]); +kernel_id = astra.data2d.create('-sino', kernel_geom, kernel); + +cfg['FilterType'] = 'rprojection' +cfg['FilterSinogramId'] = kernel_id +alg_id = astra.algorithm.create(cfg) +astra.algorithm.run(alg_id) +rec_kernel = astra.data2d.get(rec_id) +astra.algorithm.delete(alg_id) + +import pylab +pylab.figure() +pylab.imshow(P) +pylab.figure() +pylab.imshow(rec_RL) +pylab.figure() +pylab.imshow(rec_filter) +pylab.figure() +pylab.imshow(rec_kernel) +pylab.show() + +astra.data2d.delete(rec_id) +astra.data2d.delete(sinogram_id) +astra.projector.delete(proj_id) + -- cgit v1.2.3