summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
authorWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2016-10-07 16:36:28 +0200
committerWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2016-10-07 16:36:28 +0200
commitf6aa2db83dfea89f9d2cfc6fcbd3da141ee77e02 (patch)
tree6de18fc05a9ee0cf033675aadfd8514a97d9e7d1 /src
parent741675b458bf23cf437a6c56f95483c5c66af774 (diff)
parent8f2b55a66db9747419e75dae5973281a7536b934 (diff)
downloadastra-f6aa2db83dfea89f9d2cfc6fcbd3da141ee77e02.tar.gz
astra-f6aa2db83dfea89f9d2cfc6fcbd3da141ee77e02.tar.bz2
astra-f6aa2db83dfea89f9d2cfc6fcbd3da141ee77e02.tar.xz
astra-f6aa2db83dfea89f9d2cfc6fcbd3da141ee77e02.zip
Merge branch 'master' into parallel_vec
Diffstat (limited to 'src')
-rw-r--r--src/CompositeGeometryManager.cpp371
-rw-r--r--src/CudaFDKAlgorithm3D.cpp20
-rw-r--r--src/Utilities.cpp34
3 files changed, 271 insertions, 154 deletions
diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp
index 084ba8c..7c4f8e6 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,25 @@ bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeDat
return doJobs(L);
}
+
+bool CCompositeGeometryManager::doFDK(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData,
+ CFloat32ProjectionData3DMemory *pProjData, bool bShortScan)
+{
+ 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;
+
+ 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 +1252,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 +1263,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 +1275,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);
+ 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 b5ce545..c7c8ed5 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"
@@ -81,6 +82,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;
}
@@ -195,6 +207,7 @@ void CCudaFDKAlgorithm3D::run(int _iNrIterations)
ASTRA_ASSERT(pReconMem);
+#if 0
bool ok = true;
ok = astraCudaFDK(pReconMem->getData(), pSinoMem->getDataConst(),
@@ -202,6 +215,13 @@ void CCudaFDKAlgorithm3D::run(int _iNrIterations)
m_bShortScan, m_iGPUIndex, m_iVoxelSuperSampling);
ASTRA_ASSERT(ok);
+#endif
+
+ CCompositeGeometryManager cgm;
+
+ cgm.doFDK(m_pProjector, pReconMem, pSinoMem, m_bShortScan);
+
+
}
//----------------------------------------------------------------------------------------
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);
+}
}