Skip to content

Commit

Permalink
Merge pull request #160 from wjp/FBP_filters
Browse files Browse the repository at this point in the history 
Custom filter support for CPU FBP
  • Loading branch information
@wjp
wjp committed Jul 18, 2018
2 parents 0d06afc + 4d741fc commit 93612c3
Show file tree
Hide file tree
Showing 19 changed files with 1,050 additions and 657 deletions.
3 changes: 2 additions & 1 deletion astra_vc14.vcxproj
View file
Expand Up @@ -509,6 +509,7 @@
<ClCompile Include="src\FanFlatProjectionGeometry2D.cpp" />
<ClCompile Include="src\FanFlatVecProjectionGeometry2D.cpp" />
<ClCompile Include="src\FilteredBackProjectionAlgorithm.cpp" />
<ClCompile Include="src\Filters.cpp" />
<ClCompile Include="src\Float32Data.cpp" />
<ClCompile Include="src\Float32Data2D.cpp" />
<ClCompile Include="src\Float32Data3D.cpp" />
Expand Down Expand Up @@ -596,6 +597,7 @@
<ClInclude Include="include\astra\FanFlatProjectionGeometry2D.h" />
<ClInclude Include="include\astra\FanFlatVecProjectionGeometry2D.h" />
<ClInclude Include="include\astra\FilteredBackProjectionAlgorithm.h" />
<ClInclude Include="include\astra\Filters.h" />
<ClInclude Include="include\astra\Float32Data.h" />
<ClInclude Include="include\astra\Float32Data2D.h" />
<ClInclude Include="include\astra\Float32Data3D.h" />
Expand Down Expand Up @@ -656,7 +658,6 @@
<ClInclude Include="include\astra\cuda\2d\fan_bp.h" />
<ClInclude Include="include\astra\cuda\2d\fan_fp.h" />
<ClInclude Include="include\astra\cuda\2d\fbp.h" />
<ClInclude Include="include\astra\cuda\2d\fbp_filters.h" />
<ClInclude Include="include\astra\cuda\2d\fft.h" />
<ClInclude Include="include\astra\cuda\2d\par_bp.h" />
<ClInclude Include="include\astra\cuda\2d\par_fp.h" />
Expand Down
9 changes: 6 additions & 3 deletions astra_vc14.vcxproj.filters
View file
Expand Up @@ -168,6 +168,9 @@
<ClCompile Include="src\Config.cpp">
<Filter>Global &amp; Other\source</Filter>
</ClCompile>
<ClCompile Include="src\Filters.cpp">
<Filter>Global &amp; Other\source</Filter>
</ClCompile>
<ClCompile Include="src\Fourier.cpp">
<Filter>Global &amp; Other\source</Filter>
</ClCompile>
Expand Down Expand Up @@ -434,6 +437,9 @@
<ClInclude Include="include\astra\Config.h">
<Filter>Global &amp; Other\headers</Filter>
</ClInclude>
<ClInclude Include="include\astra\Filters.h">
<Filter>Global &amp; Other\headers</Filter>
</ClInclude>
<ClInclude Include="include\astra\Fourier.h">
<Filter>Global &amp; Other\headers</Filter>
</ClInclude>
Expand Down Expand Up @@ -638,9 +644,6 @@
<ClInclude Include="include\astra\cuda\2d\fan_fp.h">
<Filter>CUDA\cuda headers</Filter>
</ClInclude>
<ClInclude Include="include\astra\cuda\2d\fbp_filters.h">
<Filter>CUDA\cuda headers</Filter>
</ClInclude>
<ClInclude Include="include\astra\cuda\2d\fbp.h">
<Filter>CUDA\cuda headers</Filter>
</ClInclude>
Expand Down
1 change: 1 addition & 0 deletions build/linux/Makefile.in
View file
Expand Up @@ -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 \
Expand Down
3 changes: 2 additions & 1 deletion build/msvc/gen.py
View file
Expand Up @@ -214,6 +214,7 @@ def create_mex_project(name, uuid14, uuid11, uuid09):
"src\\AstraObjectManager.cpp",
"src\\CompositeGeometryManager.cpp",
"src\\Config.cpp",
"src\\Filters.cpp",
"src\\Fourier.cpp",
"src\\Globals.cpp",
"src\\Logging.cpp",
Expand Down Expand Up @@ -292,7 +293,6 @@ def create_mex_project(name, uuid14, uuid11, uuid09):
"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",
Expand Down Expand Up @@ -354,6 +354,7 @@ def create_mex_project(name, uuid14, uuid11, uuid09):
"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",
Expand Down
55 changes: 23 additions & 32 deletions cuda/2d/fbp.cu
View file
Expand Up @@ -35,27 +35,18 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#include "astra/cuda/3d/fdk.h"

#include "astra/Logging.h"
#include "astra/Filters.h"

#include <cuda.h>

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 iFreqBinCount = calcFFTFourierSize(iFFTRealDetCount);
int iFFTRealDetCount = astra::calcNextPowerOfTwo(2 * _iDetectorCount);
int iFreqBinCount = astra::calcFFTFourierSize(iFFTRealDetCount);

// CHECKME: Matlab makes this at least 64. Do we also need to?
return iFreqBinCount;
Expand Down Expand Up @@ -88,27 +79,27 @@ 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)
{
freeComplexOnDevice((cufftComplex*)D_filter);
D_filter = 0;
}

if (_eFilter == astra::FILTER_NONE)
if (_cfg.m_eType == astra::FILTER_NONE)
return true; // leave D_filter set to 0


int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets);
int iFreqBinCount = calcFFTFourierSize(iFFTRealDetCount);
int iFFTRealDetCount = astra::calcNextPowerOfTwo(2 * dims.iProjDets);
int iFreqBinCount = astra::calcFFTFourierSize(iFFTRealDetCount);

cufftComplex * pHostFilter = new cufftComplex[dims.iProjAngles * iFreqBinCount];
memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iFreqBinCount);

allocateComplexOnDevice(dims.iProjAngles, iFreqBinCount, (cufftComplex**)&D_filter);

switch(_eFilter)
switch(_cfg.m_eType)
{
case astra::FILTER_NONE:
// handled above
Expand All @@ -130,19 +121,19 @@ 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(_cfg, dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount);
uploadComplexArrayToDevice(dims.iProjAngles, iFreqBinCount, pHostFilter, (cufftComplex*)D_filter);

break;
}
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++)
{
Expand All @@ -156,13 +147,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;
Expand All @@ -178,16 +169,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++)
{
Expand All @@ -213,19 +204,19 @@ 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;

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;
}
}
Expand Down Expand Up @@ -310,8 +301,8 @@ bool FBP::iterate(unsigned int iterations)

if (D_filter) {

int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets);
int iFFTFourDetCount = calcFFTFourierSize(iFFTRealDetCount);
int iFFTRealDetCount = astra::calcNextPowerOfTwo(2 * dims.iProjDets);
int iFFTFourDetCount = astra::calcFFTFourierSize(iFFTRealDetCount);

cufftComplex * pDevComplexSinogram = NULL;

Expand Down

0 comments on commit 93612c3

Please sign in to comment.