/* ----------------------------------------------------------------------- Copyright: 2010-2021, imec Vision Lab, University of Antwerp 2014-2021, 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 "astra/Config.h" #include "astra/Float32ProjectionData2D.h" #include "astra/AstraObjectManager.h" #include #include namespace astra { float *genFilter(const SFilterConfig &_cfg, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount) { 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 = 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] = PI * 2.0f * fRelIndex; } switch(_cfg.m_eType) { 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 / _cfg.m_fD) / (pfW[iDetectorIndex] / 2.0f / _cfg.m_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 / _cfg.m_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] / _cfg.m_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] / _cfg.m_fD)) / 2.0f; } break; } case FILTER_TUKEY: { 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; 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 = 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 = 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 = _cfg.m_fParameter; if(_cfg.m_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 * PI * fSmallN / fNMinusOne; float fSecondTerm = fA2 * cosf(fCosInput); float fStoredValue = fA0 - fFirstTerm - fSecondTerm; pfFilt[iDetectorIndex] *= fStoredValue; } break; } case FILTER_BLACKMAN: { 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; float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1); for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { float fSmallN = (float)iDetectorIndex; float fCosInput1 = 2.0f * PI * 0.5f * fSmallN / fNMinusOne; float fCosInput2 = 4.0f * 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 = 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 = 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 = 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 = 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 = _cfg.m_fParameter; if(_cfg.m_fParameter < 0.0f) fAlpha = 3.0f; float fPiTimesAlpha = 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 = PI * _cfg.m_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 std::string &_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.c_str(), i->m_name)) return i->m_type; ASTRA_ERROR("Failed to convert \"%s\" into a filter.",_filterType.c_str()); return FILTER_ERROR; } SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg) { ConfigStackCheck CC("getFilterConfig", _alg, _cfg); SFilterConfig c; XMLNode node; // filter type const char *nodeName = "FilterType"; node = _cfg.self.getSingleNode(nodeName); if (_cfg.self.hasOption(nodeName)) { c.m_eType = convertStringToFilter(_cfg.self.getOption(nodeName)); CC.markOptionParsed(nodeName); } else if (node) { // Fallback: check cfg.FilterType (instead of cfg.option.FilterType) c.m_eType = convertStringToFilter(node.getContent()); CC.markNodeParsed(nodeName); } else { c.m_eType = FILTER_RAMLAK; } // filter nodeName = "FilterSinogramId"; int id = -1; switch (c.m_eType) { case FILTER_PROJECTION: case FILTER_SINOGRAM: case FILTER_RPROJECTION: case FILTER_RSINOGRAM: node = _cfg.self.getSingleNode(nodeName); try { if (_cfg.self.hasOption(nodeName)) { id = _cfg.self.getOptionInt(nodeName); CC.markOptionParsed(nodeName); } else if (node) { id = node.getContentInt(); CC.markNodeParsed(nodeName); } } catch (const astra::StringUtil::bad_cast &e) { ASTRA_ERROR("%s is not a valid id", nodeName); } break; default: break; } if (id != -1) { const CFloat32ProjectionData2D * pFilterData = dynamic_cast(CData2DManager::getSingleton().get(id)); c.m_iCustomFilterWidth = pFilterData->getGeometry()->getDetectorCount(); c.m_iCustomFilterHeight = pFilterData->getGeometry()->getProjectionAngleCount(); 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; } // filter parameter nodeName = "FilterParameter"; c.m_fParameter = -1.0f; switch (c.m_eType) { case FILTER_TUKEY: case FILTER_GAUSSIAN: case FILTER_BLACKMAN: case FILTER_KAISER: try { node = _cfg.self.getSingleNode(nodeName); if (_cfg.self.hasOption(nodeName)) { c.m_fParameter = _cfg.self.getOptionNumerical(nodeName); CC.markOptionParsed(nodeName); } else if (node) { c.m_fParameter = node.getContentNumerical(); CC.markNodeParsed(nodeName); } } catch (const astra::StringUtil::bad_cast &e) { ASTRA_ERROR("%s must be numerical", nodeName); } break; default: break; } // D value nodeName = "FilterD"; c.m_fD = 1.0f; switch (c.m_eType) { case FILTER_PROJECTION: case FILTER_SINOGRAM: case FILTER_RPROJECTION: case FILTER_RSINOGRAM: break; case FILTER_NONE: case FILTER_ERROR: break; default: try { node = _cfg.self.getSingleNode(nodeName); if (_cfg.self.hasOption(nodeName)) { c.m_fD = _cfg.self.getOptionNumerical(nodeName); CC.markOptionParsed(nodeName); } else if (node) { c.m_fD = node.getContentNumerical(); CC.markNodeParsed(nodeName); } } catch (const astra::StringUtil::bad_cast &e) { ASTRA_ERROR("%s must be numerical", nodeName); } break; } return c; } int calcNextPowerOfTwo(int n) { int x = 1; while (x < n && x > 0) 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; } } }