diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/CompositeGeometryManager.cpp | 373 | ||||
| -rw-r--r-- | src/CudaFDKAlgorithm3D.cpp | 32 | ||||
| -rw-r--r-- | src/Utilities.cpp | 34 |
3 files changed, 279 insertions, 160 deletions
diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index 084ba8c..5879aec 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -146,6 +146,7 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div newjob.eType = j->eType; newjob.eMode = j->eMode; newjob.pProjector = j->pProjector; + newjob.FDKSettings = j->FDKSettings; CPart* input = job.pInput->reduce(outputPart.get()); @@ -196,6 +197,69 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div return true; } + +static std::pair<double, double> reduceProjectionVertical(const CVolumeGeometry3D* pVolGeom, const CProjectionGeometry3D* pProjGeom) +{ + double vmin_g, vmax_g; + + // reduce self to only cover intersection with projection of VolumePart + // (Project corners of volume, take bounding box) + + for (int i = 0; i < pProjGeom->getProjectionCount(); ++i) { + + double vol_u[8]; + double vol_v[8]; + + double pixx = pVolGeom->getPixelLengthX(); + double pixy = pVolGeom->getPixelLengthY(); + double pixz = pVolGeom->getPixelLengthZ(); + + // TODO: Is 0.5 sufficient? + double xmin = pVolGeom->getWindowMinX() - 0.5 * pixx; + double xmax = pVolGeom->getWindowMaxX() + 0.5 * pixx; + double ymin = pVolGeom->getWindowMinY() - 0.5 * pixy; + double ymax = pVolGeom->getWindowMaxY() + 0.5 * pixy; + double zmin = pVolGeom->getWindowMinZ() - 0.5 * pixz; + double zmax = pVolGeom->getWindowMaxZ() + 0.5 * pixz; + + pProjGeom->projectPoint(xmin, ymin, zmin, i, vol_u[0], vol_v[0]); + pProjGeom->projectPoint(xmin, ymin, zmax, i, vol_u[1], vol_v[1]); + pProjGeom->projectPoint(xmin, ymax, zmin, i, vol_u[2], vol_v[2]); + pProjGeom->projectPoint(xmin, ymax, zmax, i, vol_u[3], vol_v[3]); + pProjGeom->projectPoint(xmax, ymin, zmin, i, vol_u[4], vol_v[4]); + pProjGeom->projectPoint(xmax, ymin, zmax, i, vol_u[5], vol_v[5]); + pProjGeom->projectPoint(xmax, ymax, zmin, i, vol_u[6], vol_v[6]); + pProjGeom->projectPoint(xmax, ymax, zmax, i, vol_u[7], vol_v[7]); + + double vmin = vol_v[0]; + double vmax = vol_v[0]; + + for (int j = 1; j < 8; ++j) { + if (vol_v[j] < vmin) + vmin = vol_v[j]; + if (vol_v[j] > vmax) + vmax = vol_v[j]; + } + + if (i == 0 || vmin < vmin_g) + vmin_g = vmin; + if (i == 0 || vmax > vmax_g) + vmax_g = vmax; + } + + if (vmin_g < -1.0) + vmin_g = -1.0; + if (vmax_g > pProjGeom->getDetectorRowCount()) + vmax_g = pProjGeom->getDetectorRowCount(); + + if (vmin_g >= vmax_g) + vmin_g = vmax_g = 0.0; + + return std::pair<double, double>(vmin_g, vmax_g); + +} + + CCompositeGeometryManager::CPart::CPart(const CPart& other) { eType = other.eType; @@ -236,119 +300,135 @@ size_t CCompositeGeometryManager::CPart::getSize() } +static bool testVolumeRange(const std::pair<double, double>& fullRange, + const CVolumeGeometry3D *pVolGeom, + const CProjectionGeometry3D *pProjGeom, + int zmin, int zmax) +{ + double pixz = pVolGeom->getPixelLengthZ(); + + CVolumeGeometry3D test(pVolGeom->getGridColCount(), + pVolGeom->getGridRowCount(), + zmax - zmin, + pVolGeom->getWindowMinX(), + pVolGeom->getWindowMinY(), + pVolGeom->getWindowMinZ() + zmin * pixz, + pVolGeom->getWindowMaxX(), + pVolGeom->getWindowMaxY(), + pVolGeom->getWindowMinZ() + zmax * pixz); + + + std::pair<double, double> subRange = reduceProjectionVertical(&test, pProjGeom); + + // empty + if (subRange.first == subRange.second) + return true; + + // fully outside of fullRange + if (subRange.first >= fullRange.second || subRange.second <= fullRange.first) + return true; + + return false; +} + CCompositeGeometryManager::CPart* CCompositeGeometryManager::CVolumePart::reduce(const CPart *_other) { const CProjectionPart *other = dynamic_cast<const CProjectionPart *>(_other); assert(other); - // TODO: Is 0.5 sufficient? - double umin = -0.5; - double umax = other->pGeom->getDetectorColCount() + 0.5; - double vmin = -0.5; - double vmax = other->pGeom->getDetectorRowCount() + 0.5; - - double uu[4]; - double vv[4]; - uu[0] = umin; vv[0] = vmin; - uu[1] = umin; vv[1] = vmax; - uu[2] = umax; vv[2] = vmin; - uu[3] = umax; vv[3] = vmax; - - double pixx = pGeom->getPixelLengthX(); - double pixy = pGeom->getPixelLengthY(); - double pixz = pGeom->getPixelLengthZ(); - double xmin = pGeom->getWindowMinX() - 0.5 * pixx; - double xmax = pGeom->getWindowMaxX() + 0.5 * pixx; - double ymin = pGeom->getWindowMinY() - 0.5 * pixy; - double ymax = pGeom->getWindowMaxY() + 0.5 * pixy; - - // NB: Flipped - double zmax = pGeom->getWindowMinZ() - 2.5 * pixz; - double zmin = pGeom->getWindowMaxZ() + 2.5 * pixz; - - // TODO: This isn't as tight as it could be. - // In particular it won't detect the detector being - // missed entirely on the u side. - - for (int i = 0; i < other->pGeom->getProjectionCount(); ++i) { - for (int j = 0; j < 4; ++j) { - double px, py, pz; - - other->pGeom->backprojectPointX(i, uu[j], vv[j], xmin, py, pz); - //ASTRA_DEBUG("%f %f (%f - %f)", py, pz, ymin, ymax); - if (pz < zmin) zmin = pz; - if (pz > zmax) zmax = pz; - other->pGeom->backprojectPointX(i, uu[j], vv[j], xmax, py, pz); - //ASTRA_DEBUG("%f %f (%f - %f)", py, pz, ymin, ymax); - if (pz < zmin) zmin = pz; - if (pz > zmax) zmax = pz; - - other->pGeom->backprojectPointY(i, uu[j], vv[j], ymin, px, pz); - //ASTRA_DEBUG("%f %f (%f - %f)", px, pz, xmin, xmax); - if (pz < zmin) zmin = pz; - if (pz > zmax) zmax = pz; - other->pGeom->backprojectPointY(i, uu[j], vv[j], ymax, px, pz); - //ASTRA_DEBUG("%f %f (%f - %f)", px, pz, xmin, xmax); - if (pz < zmin) zmin = pz; - if (pz > zmax) zmax = pz; + std::pair<double, double> fullRange = reduceProjectionVertical(pGeom, other->pGeom); + + int top_slice = 0, bottom_slice = 0; + + if (fullRange.first < fullRange.second) { + + + // TOP SLICE + + int zmin = 0; + int zmax = pGeom->getGridSliceCount()-1; // (Don't try empty region) + + // Setting top slice to zmin is always valid. + + while (zmin < zmax) { + int zmid = (zmin + zmax + 1) / 2; + + bool ok = testVolumeRange(fullRange, pGeom, other->pGeom, + 0, zmid); + + ASTRA_DEBUG("binsearch min: [%d,%d], %d, %s", zmin, zmax, zmid, ok ? "ok" : "removed too much"); + + if (ok) + zmin = zmid; + else + zmax = zmid - 1; } - } - //ASTRA_DEBUG("coord extent: %f - %f", zmin, zmax); + top_slice = zmin; + + + // BOTTOM SLICE + + zmin = top_slice + 1; // (Don't try empty region) + zmax = pGeom->getGridSliceCount(); + + // Setting bottom slice to zmax is always valid - // Clip both zmin and zmax to get rid of extreme (or infinite) values - // NB: When individual pz values are +/-Inf, the sign is determined - // by ray direction and on which side of the face the ray passes. - if (zmin < pGeom->getWindowMinZ() - 2*pixz) - zmin = pGeom->getWindowMinZ() - 2*pixz; - if (zmin > pGeom->getWindowMaxZ() + 2*pixz) - zmin = pGeom->getWindowMaxZ() + 2*pixz; - if (zmax < pGeom->getWindowMinZ() - 2*pixz) - zmax = pGeom->getWindowMinZ() - 2*pixz; - if (zmax > pGeom->getWindowMaxZ() + 2*pixz) - zmax = pGeom->getWindowMaxZ() + 2*pixz; + while (zmin < zmax) { + int zmid = (zmin + zmax) / 2; - zmin = (zmin - pixz - pGeom->getWindowMinZ()) / pixz; - zmax = (zmax + pixz - pGeom->getWindowMinZ()) / pixz; + bool ok = testVolumeRange(fullRange, pGeom, other->pGeom, + zmid, pGeom->getGridSliceCount()); - int _zmin = (int)floor(zmin); - int _zmax = (int)ceil(zmax); + ASTRA_DEBUG("binsearch max: [%d,%d], %d, %s", zmin, zmax, zmid, ok ? "ok" : "removed too much"); - //ASTRA_DEBUG("index extent: %d - %d", _zmin, _zmax); + if (ok) + zmax = zmid; + else + zmin = zmid + 1; - if (_zmin < 0) - _zmin = 0; - if (_zmax > pGeom->getGridSliceCount()) - _zmax = pGeom->getGridSliceCount(); + } + + bottom_slice = zmax; - if (_zmax <= _zmin) { - _zmin = _zmax = 0; } - //ASTRA_DEBUG("adjusted extent: %d - %d", _zmin, _zmax); + + ASTRA_DEBUG("found extent: %d - %d", top_slice, bottom_slice); + + top_slice -= 1; + if (top_slice < 0) + top_slice = 0; + bottom_slice += 1; + if (bottom_slice >= pGeom->getGridSliceCount()) + bottom_slice = pGeom->getGridSliceCount(); + + ASTRA_DEBUG("adjusted extent: %d - %d", top_slice, bottom_slice); + + double pixz = pGeom->getPixelLengthZ(); CVolumePart *sub = new CVolumePart(); sub->subX = this->subX; sub->subY = this->subY; - sub->subZ = this->subZ + _zmin; + sub->subZ = this->subZ + top_slice; sub->pData = pData; - if (_zmin == _zmax) { + if (top_slice == bottom_slice) { sub->pGeom = 0; } else { sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(), pGeom->getGridRowCount(), - _zmax - _zmin, + bottom_slice - top_slice, pGeom->getWindowMinX(), pGeom->getWindowMinY(), - pGeom->getWindowMinZ() + _zmin * pixz, + pGeom->getWindowMinZ() + top_slice * pixz, pGeom->getWindowMaxX(), pGeom->getWindowMaxY(), - pGeom->getWindowMinZ() + _zmax * pixz); + pGeom->getWindowMinZ() + bottom_slice * pixz); } - ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d", this->subZ, this->subZ + pGeom->getGridSliceCount(), this->subZ + _zmin, this->subZ + _zmax); + ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d ( %f - %f )", this->subZ, this->subZ + pGeom->getGridSliceCount(), this->subZ + top_slice, this->subZ + bottom_slice, pGeom->getWindowMinZ() + top_slice * pixz, pGeom->getWindowMinZ() + bottom_slice * pixz); return sub; } @@ -583,7 +663,9 @@ void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::T size_t m = std::min(maxSize / sliceSize, maxDim); size_t blockSize = computeLinearSplit(m, div, sliceCount); - int rem = sliceCount % blockSize; + int rem = blockSize - (sliceCount % blockSize); + if (rem == blockSize) + rem = 0; ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); @@ -630,7 +712,9 @@ void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::T size_t m = std::min(maxSize / sliceSize, maxDim); size_t blockSize = computeLinearSplit(m, div, sliceCount); - int rem = sliceCount % blockSize; + int rem = blockSize - (sliceCount % blockSize); + if (rem == blockSize) + rem = 0; ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); @@ -677,7 +761,9 @@ void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::T size_t m = std::min(maxSize / sliceSize, maxDim); size_t blockSize = computeLinearSplit(m, div, sliceCount); - int rem = sliceCount % blockSize; + int rem = blockSize - (sliceCount % blockSize); + if (rem == blockSize) + rem = 0; ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); @@ -742,62 +828,16 @@ void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, s } + CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::reduce(const CPart *_other) { const CVolumePart *other = dynamic_cast<const CVolumePart *>(_other); assert(other); - double vmin_g, vmax_g; - - // reduce self to only cover intersection with projection of VolumePart - // (Project corners of volume, take bounding box) - - for (int i = 0; i < pGeom->getProjectionCount(); ++i) { - - double vol_u[8]; - double vol_v[8]; - - double pixx = other->pGeom->getPixelLengthX(); - double pixy = other->pGeom->getPixelLengthY(); - double pixz = other->pGeom->getPixelLengthZ(); - - // TODO: Is 0.5 sufficient? - double xmin = other->pGeom->getWindowMinX() - 0.5 * pixx; - double xmax = other->pGeom->getWindowMaxX() + 0.5 * pixx; - double ymin = other->pGeom->getWindowMinY() - 0.5 * pixy; - double ymax = other->pGeom->getWindowMaxY() + 0.5 * pixy; - double zmin = other->pGeom->getWindowMinZ() - 0.5 * pixz; - double zmax = other->pGeom->getWindowMaxZ() + 0.5 * pixz; - - pGeom->projectPoint(xmin, ymin, zmin, i, vol_u[0], vol_v[0]); - pGeom->projectPoint(xmin, ymin, zmax, i, vol_u[1], vol_v[1]); - pGeom->projectPoint(xmin, ymax, zmin, i, vol_u[2], vol_v[2]); - pGeom->projectPoint(xmin, ymax, zmax, i, vol_u[3], vol_v[3]); - pGeom->projectPoint(xmax, ymin, zmin, i, vol_u[4], vol_v[4]); - pGeom->projectPoint(xmax, ymin, zmax, i, vol_u[5], vol_v[5]); - pGeom->projectPoint(xmax, ymax, zmin, i, vol_u[6], vol_v[6]); - pGeom->projectPoint(xmax, ymax, zmax, i, vol_u[7], vol_v[7]); - - double vmin = vol_v[0]; - double vmax = vol_v[0]; - - for (int j = 1; j < 8; ++j) { - if (vol_v[j] < vmin) - vmin = vol_v[j]; - if (vol_v[j] > vmax) - vmax = vol_v[j]; - } - - if (i == 0 || vmin < vmin_g) - vmin_g = vmin; - if (i == 0 || vmax > vmax_g) - vmax_g = vmax; - } - - // fprintf(stderr, "v extent: %f %f\n", vmin_g, vmax_g); - - int _vmin = (int)floor(vmin_g - 1.0f); - int _vmax = (int)ceil(vmax_g + 1.0f); + std::pair<double, double> r = reduceProjectionVertical(other->pGeom, pGeom); + // fprintf(stderr, "v extent: %f %f\n", r.first, r.second); + int _vmin = (int)floor(r.first - 1.0); + int _vmax = (int)ceil(r.second + 1.0); if (_vmin < 0) _vmin = 0; if (_vmax > pGeom->getDetectorRowCount()) @@ -837,7 +877,11 @@ void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManage size_t m = std::min(maxSize / sliceSize, maxDim); size_t blockSize = computeLinearSplit(m, div, sliceCount); - int rem = sliceCount % blockSize; + int rem = blockSize - (sliceCount % blockSize); + if (rem == blockSize) + rem = 0; + + ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); for (int x = -(rem / 2); x < sliceCount; x += blockSize) { int newsubX = x; @@ -879,7 +923,11 @@ void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManage size_t m = std::min(maxSize / sliceSize, maxDim); size_t blockSize = computeLinearSplit(m, div, sliceCount); - int rem = sliceCount % blockSize; + int rem = blockSize - (sliceCount % blockSize); + if (rem == blockSize) + rem = 0; + + ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); for (int z = -(rem / 2); z < sliceCount; z += blockSize) { int newsubZ = z; @@ -992,6 +1040,27 @@ bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeDat return doJobs(L); } + +bool CCompositeGeometryManager::doFDK(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, + CFloat32ProjectionData3DMemory *pProjData, bool bShortScan, + const float *pfFilter) +{ + if (!dynamic_cast<CConeProjectionGeometry3D*>(pProjData->getGeometry())) { + ASTRA_ERROR("CCompositeGeometryManager::doFDK: cone geometry required"); + return false; + } + + SJob job = createJobBP(pProjector, pVolData, pProjData); + job.eType = SJob::JOB_FDK; + job.FDKSettings.bShortScan = bShortScan; + job.FDKSettings.pfFilter = pfFilter; + + TJobList L; + L.push_back(job); + + return doJobs(L); +} + bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData) { ASTRA_DEBUG("CCompositeGeometryManager::doFP, multi-volume"); @@ -1185,7 +1254,9 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter ok = astraCUDA3d::copyToGPUMemory(src, inputMem, srcdims); if (!ok) ASTRA_ERROR("Error copying input data to GPU"); - if (j.eType == CCompositeGeometryManager::SJob::JOB_FP) { + switch (j.eType) { + case CCompositeGeometryManager::SJob::JOB_FP: + { assert(dynamic_cast<CCompositeGeometryManager::CVolumePart*>(j.pInput.get())); assert(dynamic_cast<CCompositeGeometryManager::CProjectionPart*>(j.pOutput.get())); @@ -1194,7 +1265,10 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter ok = astraCUDA3d::FP(((CCompositeGeometryManager::CProjectionPart*)j.pOutput.get())->pGeom, outputMem, ((CCompositeGeometryManager::CVolumePart*)j.pInput.get())->pGeom, inputMem, detectorSuperSampling, projKernel); if (!ok) ASTRA_ERROR("Error performing sub-FP"); ASTRA_DEBUG("CCompositeGeometryManager::doJobs: FP done"); - } else if (j.eType == CCompositeGeometryManager::SJob::JOB_BP) { + } + break; + case CCompositeGeometryManager::SJob::JOB_BP: + { assert(dynamic_cast<CCompositeGeometryManager::CVolumePart*>(j.pOutput.get())); assert(dynamic_cast<CCompositeGeometryManager::CProjectionPart*>(j.pInput.get())); @@ -1203,7 +1277,26 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter ok = astraCUDA3d::BP(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, inputMem, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, outputMem, voxelSuperSampling); if (!ok) ASTRA_ERROR("Error performing sub-BP"); ASTRA_DEBUG("CCompositeGeometryManager::doJobs: BP done"); - } else { + } + break; + case CCompositeGeometryManager::SJob::JOB_FDK: + { + assert(dynamic_cast<CCompositeGeometryManager::CVolumePart*>(j.pOutput.get())); + assert(dynamic_cast<CCompositeGeometryManager::CProjectionPart*>(j.pInput.get())); + + if (srcdims.subx || srcdims.suby) { + ASTRA_ERROR("CCompositeGeometryManager::doJobs: data too large for FDK"); + ok = false; + } else { + ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing FDK"); + + ok = astraCUDA3d::FDK(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, inputMem, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, outputMem, j.FDKSettings.bShortScan, j.FDKSettings.pfFilter); + if (!ok) ASTRA_ERROR("Error performing sub-FDK"); + ASTRA_DEBUG("CCompositeGeometryManager::doJobs: FDK done"); + } + } + break; + default: assert(false); } diff --git a/src/CudaFDKAlgorithm3D.cpp b/src/CudaFDKAlgorithm3D.cpp index bc15a3d..15abed4 100644 --- a/src/CudaFDKAlgorithm3D.cpp +++ b/src/CudaFDKAlgorithm3D.cpp @@ -32,6 +32,7 @@ $Id$ #include "astra/CudaProjector3D.h" #include "astra/ConeProjectionGeometry3D.h" +#include "astra/CompositeGeometryManager.h" #include "astra/Logging.h" @@ -84,6 +85,17 @@ bool CCudaFDKAlgorithm3D::_check() const CProjectionGeometry3D* projgeom = m_pSinogram->getGeometry(); ASTRA_CONFIG_CHECK(dynamic_cast<const CConeProjectionGeometry3D*>(projgeom), "CUDA_FDK", "Error setting FDK geometry"); + + const CVolumeGeometry3D* volgeom = m_pReconstruction->getGeometry(); + bool cube = true; + if (abs(volgeom->getPixelLengthX() / volgeom->getPixelLengthY() - 1.0) > 0.00001) + cube = false; + if (abs(volgeom->getPixelLengthX() / volgeom->getPixelLengthZ() - 1.0) > 0.00001) + cube = false; + ASTRA_CONFIG_CHECK(cube, "CUDA_FDK", "Voxels must be cubes for FDK"); + + + return true; } @@ -219,20 +231,28 @@ void CCudaFDKAlgorithm3D::run(int _iNrIterations) CFloat32VolumeData3DMemory* pReconMem = dynamic_cast<CFloat32VolumeData3DMemory*>(m_pReconstruction); ASTRA_ASSERT(pReconMem); - - bool ok = true; - const float *filter = NULL; - if(m_iFilterDataId!=-1){ - const CFloat32ProjectionData2D * pFilterData = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(m_iFilterDataId)); - filter = pFilterData->getDataConst(); + if (m_iFilterDataId != -1) { + const CFloat32ProjectionData2D *pFilterData = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(m_iFilterDataId)); + if (pFilterData) + filter = pFilterData->getDataConst(); } +#if 0 + bool ok = true; + ok = astraCudaFDK(pReconMem->getData(), pSinoMem->getDataConst(), &volgeom, conegeom, m_bShortScan, m_iGPUIndex, m_iVoxelSuperSampling, filter); ASTRA_ASSERT(ok); +#endif + + CCompositeGeometryManager cgm; + + cgm.doFDK(m_pProjector, pReconMem, pSinoMem, m_bShortScan, filter); + + } //---------------------------------------------------------------------------------------- diff --git a/src/Utilities.cpp b/src/Utilities.cpp index c9740bf..8b0ca94 100644 --- a/src/Utilities.cpp +++ b/src/Utilities.cpp @@ -28,10 +28,6 @@ $Id$ #include "astra/Utilities.h" -#include <boost/algorithm/string.hpp> -#include <boost/algorithm/string/split.hpp> -#include <boost/algorithm/string/classification.hpp> - #include <sstream> #include <locale> #include <iomanip> @@ -84,18 +80,16 @@ std::vector<double> stringToDoubleVector(const std::string &s) template<typename T> std::vector<T> stringToVector(const std::string& s) { - // split - std::vector<std::string> items; - boost::split(items, s, boost::is_any_of(",;")); - - // init list std::vector<T> out; - out.resize(items.size()); + size_t current = 0; + size_t next; + do { + next = s.find_first_of(",;", current); + std::string t = s.substr(current, next - current); + out.push_back(stringTo<T>(t)); + current = next + 1; + } while (next != std::string::npos); - // loop elements - for (unsigned int i = 0; i < items.size(); i++) { - out[i] = stringTo<T>(items[i]); - } return out; } @@ -120,6 +114,18 @@ std::string doubleToString(double f) template<> std::string toString(float f) { return floatToString(f); } template<> std::string toString(double f) { return doubleToString(f); } +void splitString(std::vector<std::string> &items, const std::string& s, + const char *delim) +{ + items.clear(); + size_t current = 0; + size_t next; + do { + next = s.find_first_of(",;", current); + items.push_back(s.substr(current, next - current)); + current = next + 1; + } while (next != std::string::npos); +} } |