diff options
Diffstat (limited to 'src')
27 files changed, 2259 insertions, 567 deletions
diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp new file mode 100644 index 0000000..9be4797 --- /dev/null +++ b/src/CompositeGeometryManager.cpp @@ -0,0 +1,993 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp + 2014-2015, CWI, Amsterdam + +Contact: astra@uantwerpen.be +Website: http://sf.net/projects/astra-toolbox + +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 <http://www.gnu.org/licenses/>. + +----------------------------------------------------------------------- +*/ + +#include "astra/CompositeGeometryManager.h" + +#ifdef ASTRA_CUDA + +#include "astra/GeometryUtil3D.h" +#include "astra/VolumeGeometry3D.h" +#include "astra/ConeProjectionGeometry3D.h" +#include "astra/ConeVecProjectionGeometry3D.h" +#include "astra/ParallelProjectionGeometry3D.h" +#include "astra/ParallelVecProjectionGeometry3D.h" +#include "astra/Projector3D.h" +#include "astra/CudaProjector3D.h" +#include "astra/Float32ProjectionData3DMemory.h" +#include "astra/Float32VolumeData3DMemory.h" +#include "astra/Logging.h" + +#include "../cuda/3d/mem3d.h" + +#include <cstring> + +namespace astra { + +// JOB: +// +// VolumePart +// ProjectionPart +// FP-or-BP +// SET-or-ADD + + +// Running a set of jobs: +// +// [ Assume OUTPUT Parts in a single JobSet don't alias?? ] +// Group jobs by output Part +// One thread per group? + +// Automatically split parts if too large +// Performance model for odd-sized tasks? +// Automatically split parts if not enough tasks to fill available GPUs + + +// Splitting: +// Constraints: +// number of sub-parts divisible by N +// max size of sub-parts + +// For splitting on both input and output side: +// How to divide up memory? (Optimization problem; compute/benchmark) +// (First approach: 0.5/0.5) + + + +bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div, TJobSet &split) +{ + split.clear(); + + for (TJobSet::const_iterator i = jobs.begin(); i != jobs.end(); ++i) + { + CPart* pOutput = i->first; + const TJobList &L = i->second; + + // 1. Split output part + // 2. Per sub-part: + // a. reduce input part + // b. split input part + // c. create jobs for new (input,output) subparts + + TPartList splitOutput = pOutput->split(maxSize/3, div); + + for (TJobList::const_iterator j = L.begin(); j != L.end(); ++j) + { + const SJob &job = *j; + + for (TPartList::iterator i_out = splitOutput.begin(); + i_out != splitOutput.end(); ++i_out) + { + boost::shared_ptr<CPart> outputPart = *i_out; + + SJob newjob; + newjob.pOutput = outputPart; + newjob.eType = j->eType; + newjob.eMode = j->eMode; + newjob.pProjector = j->pProjector; + + CPart* input = job.pInput->reduce(outputPart.get()); + + if (input->getSize() == 0) { + ASTRA_DEBUG("Empty input"); + newjob.eType = SJob::JOB_NOP; + split[outputPart.get()].push_back(newjob); + continue; + } + + size_t remainingSize = ( maxSize - outputPart->getSize() ) / 2; + + TPartList splitInput = input->split(remainingSize, 1); + delete input; + ASTRA_DEBUG("Input split into %d parts", splitInput.size()); + + for (TPartList::iterator i_in = splitInput.begin(); + i_in != splitInput.end(); ++i_in) + { + newjob.pInput = *i_in; + + split[outputPart.get()].push_back(newjob); + + // Second and later (input) parts should always be added to + // output of first (input) part. + newjob.eMode = SJob::MODE_ADD; + } + + + } + + } + } + + return true; +} + +CCompositeGeometryManager::CPart::CPart(const CPart& other) +{ + eType = other.eType; + pData = other.pData; + subX = other.subX; + subY = other.subY; + subZ = other.subZ; +} + +CCompositeGeometryManager::CVolumePart::CVolumePart(const CVolumePart& other) + : CPart(other) +{ + pGeom = other.pGeom->clone(); +} + +CCompositeGeometryManager::CVolumePart::~CVolumePart() +{ + delete pGeom; +} + +void CCompositeGeometryManager::CVolumePart::getDims(size_t &x, size_t &y, size_t &z) +{ + if (!pGeom) { + x = y = z = 0; + return; + } + + x = pGeom->getGridColCount(); + y = pGeom->getGridRowCount(); + z = pGeom->getGridSliceCount(); +} + +size_t CCompositeGeometryManager::CPart::getSize() +{ + size_t x, y, z; + getDims(x, y, z); + return x * y * z; +} + + + +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; + } + } + + //ASTRA_DEBUG("coord extent: %f - %f", zmin, zmax); + + zmin = (zmin - pixz - pGeom->getWindowMinZ()) / pixz; + zmax = (zmax + pixz - pGeom->getWindowMinZ()) / pixz; + + int _zmin = (int)floor(zmin); + int _zmax = (int)ceil(zmax); + + //ASTRA_DEBUG("index extent: %d - %d", _zmin, _zmax); + + if (_zmin < 0) + _zmin = 0; + if (_zmax > pGeom->getGridSliceCount()) + _zmax = pGeom->getGridSliceCount(); + + if (_zmax <= _zmin) { + _zmin = _zmax = 0; + } + //ASTRA_DEBUG("adjusted extent: %d - %d", _zmin, _zmax); + + CVolumePart *sub = new CVolumePart(); + sub->subX = this->subX; + sub->subY = this->subY; + sub->subZ = this->subZ + _zmin; + sub->pData = pData; + + if (_zmin == _zmax) { + sub->pGeom = 0; + } else { + sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(), + pGeom->getGridRowCount(), + _zmax - _zmin, + pGeom->getWindowMinX(), + pGeom->getWindowMinY(), + pGeom->getWindowMinZ() + _zmin * pixz, + pGeom->getWindowMaxX(), + pGeom->getWindowMaxY(), + pGeom->getWindowMinZ() + _zmax * pixz); + } + + ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d", this->subZ, this->subZ + pGeom->getGridSliceCount(), this->subZ + _zmin, this->subZ + _zmax); + + return sub; +} + + + +static size_t ceildiv(size_t a, size_t b) { + return (a + b - 1) / b; +} + +static size_t computeVerticalSplit(size_t maxBlock, int div, size_t sliceCount) +{ + size_t blockSize = maxBlock; + size_t blockCount = ceildiv(sliceCount, blockSize); + + // Increase number of blocks to be divisible by div + size_t divCount = div * ceildiv(blockCount, div); + + // If divCount is above sqrt(number of slices), then + // we can't guarantee divisibility by div, but let's try anyway + if (ceildiv(sliceCount, ceildiv(sliceCount, divCount)) % div == 0) { + blockCount = divCount; + } else { + // If divisibility isn't achievable, we may want to optimize + // differently. + // TODO: Figure out how to model and optimize this. + } + + // Final adjustment to make blocks more evenly sized + // (This can't make the blocks larger) + blockSize = ceildiv(sliceCount, blockCount); + + ASTRA_DEBUG("%ld %ld -> %ld * %ld", sliceCount, maxBlock, blockCount, blockSize); + + assert(blockSize <= maxBlock); + assert((divCount * divCount > sliceCount) || (blockCount % div) == 0); + + return blockSize; +} + +template<class V, class P> +static V* getProjectionVectors(const P* geom); + +template<> +SConeProjection* getProjectionVectors(const CConeProjectionGeometry3D* pProjGeom) +{ + return genConeProjections(pProjGeom->getProjectionCount(), + pProjGeom->getDetectorColCount(), + pProjGeom->getDetectorRowCount(), + pProjGeom->getOriginSourceDistance(), + pProjGeom->getOriginDetectorDistance(), + pProjGeom->getDetectorSpacingX(), + pProjGeom->getDetectorSpacingY(), + pProjGeom->getProjectionAngles()); +} + +template<> +SConeProjection* getProjectionVectors(const CConeVecProjectionGeometry3D* pProjGeom) +{ + int nth = pProjGeom->getProjectionCount(); + + SConeProjection* pProjs = new SConeProjection[nth]; + for (int i = 0; i < nth; ++i) + pProjs[i] = pProjGeom->getProjectionVectors()[i]; + + return pProjs; +} + +template<> +SPar3DProjection* getProjectionVectors(const CParallelProjectionGeometry3D* pProjGeom) +{ + return genPar3DProjections(pProjGeom->getProjectionCount(), + pProjGeom->getDetectorColCount(), + pProjGeom->getDetectorRowCount(), + pProjGeom->getDetectorSpacingX(), + pProjGeom->getDetectorSpacingY(), + pProjGeom->getProjectionAngles()); +} + +template<> +SPar3DProjection* getProjectionVectors(const CParallelVecProjectionGeometry3D* pProjGeom) +{ + int nth = pProjGeom->getProjectionCount(); + + SPar3DProjection* pProjs = new SPar3DProjection[nth]; + for (int i = 0; i < nth; ++i) + pProjs[i] = pProjGeom->getProjectionVectors()[i]; + + return pProjs; +} + + +template<class V> +static void translateProjectionVectors(V* pProjs, int count, double dv) +{ + for (int i = 0; i < count; ++i) { + pProjs[i].fDetSX += dv * pProjs[i].fDetVX; + pProjs[i].fDetSY += dv * pProjs[i].fDetVY; + pProjs[i].fDetSZ += dv * pProjs[i].fDetVZ; + } +} + + + +static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry3D* pProjGeom, int v, int size) +{ + // First convert to vectors, then translate, then convert into new object + + const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(pProjGeom); + const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(pProjGeom); + const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(pProjGeom); + const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(pProjGeom); + + if (conegeom || conevec3dgeom) { + SConeProjection* pConeProjs; + if (conegeom) { + pConeProjs = getProjectionVectors<SConeProjection>(conegeom); + } else { + pConeProjs = getProjectionVectors<SConeProjection>(conevec3dgeom); + } + + translateProjectionVectors(pConeProjs, pProjGeom->getProjectionCount(), v); + + CProjectionGeometry3D* ret = new CConeVecProjectionGeometry3D(pProjGeom->getProjectionCount(), + size, + pProjGeom->getDetectorColCount(), + pConeProjs); + + + delete[] pConeProjs; + return ret; + } else { + assert(par3dgeom || parvec3dgeom); + SPar3DProjection* pParProjs; + if (par3dgeom) { + pParProjs = getProjectionVectors<SPar3DProjection>(par3dgeom); + } else { + pParProjs = getProjectionVectors<SPar3DProjection>(parvec3dgeom); + } + + translateProjectionVectors(pParProjs, pProjGeom->getProjectionCount(), v); + + CProjectionGeometry3D* ret = new CParallelVecProjectionGeometry3D(pProjGeom->getProjectionCount(), + size, + pProjGeom->getDetectorColCount(), + pParProjs); + + delete[] pParProjs; + return ret; + } + +} + + + +// split self into sub-parts: +// - each no bigger than maxSize +// - number of sub-parts is divisible by div +// - maybe all approximately the same size? +CCompositeGeometryManager::TPartList CCompositeGeometryManager::CVolumePart::split(size_t maxSize, int div) +{ + TPartList ret; + + if (true) { + // Split in vertical direction only at first, until we figure out + // a model for splitting in other directions + + size_t sliceSize = ((size_t) pGeom->getGridColCount()) * pGeom->getGridRowCount(); + int sliceCount = pGeom->getGridSliceCount(); + size_t blockSize = computeVerticalSplit(maxSize / sliceSize, div, sliceCount); + + int rem = sliceCount % blockSize; + + ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); + + for (int z = -(rem / 2); z < sliceCount; z += blockSize) { + int newsubZ = z; + if (newsubZ < 0) newsubZ = 0; + int endZ = z + blockSize; + if (endZ > sliceCount) endZ = sliceCount; + int size = endZ - newsubZ; + + CVolumePart *sub = new CVolumePart(); + sub->subX = this->subX; + sub->subY = this->subY; + sub->subZ = this->subZ + newsubZ; + + ASTRA_DEBUG("VolumePart split %d %d %d -> %p", sub->subX, sub->subY, sub->subZ, (void*)sub); + + double shift = pGeom->getPixelLengthZ() * newsubZ; + + sub->pData = pData; + sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(), + pGeom->getGridRowCount(), + size, + pGeom->getWindowMinX(), + pGeom->getWindowMinY(), + pGeom->getWindowMinZ() + shift, + pGeom->getWindowMaxX(), + pGeom->getWindowMaxY(), + pGeom->getWindowMinZ() + shift + size * pGeom->getPixelLengthZ()); + + ret.push_back(boost::shared_ptr<CPart>(sub)); + } + } + + return ret; +} + +CCompositeGeometryManager::CVolumePart* CCompositeGeometryManager::CVolumePart::clone() const +{ + return new CVolumePart(*this); +} + +CCompositeGeometryManager::CProjectionPart::CProjectionPart(const CProjectionPart& other) + : CPart(other) +{ + pGeom = other.pGeom->clone(); +} + +CCompositeGeometryManager::CProjectionPart::~CProjectionPart() +{ + delete pGeom; +} + +void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, size_t &z) +{ + if (!pGeom) { + x = y = z = 0; + return; + } + + x = pGeom->getDetectorColCount(); + y = pGeom->getProjectionCount(); + z = pGeom->getDetectorRowCount(); +} + + +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); + if (_vmin < 0) + _vmin = 0; + if (_vmax > pGeom->getDetectorRowCount()) + _vmax = pGeom->getDetectorRowCount(); + + if (_vmin >= _vmax) { + _vmin = _vmax = 0; + } + + CProjectionPart *sub = new CProjectionPart(); + sub->subX = this->subX; + sub->subY = this->subY; + sub->subZ = this->subZ + _vmin; + + sub->pData = pData; + + if (_vmin == _vmax) { + sub->pGeom = 0; + } else { + sub->pGeom = getSubProjectionGeometry(pGeom, _vmin, _vmax - _vmin); + } + + ASTRA_DEBUG("Reduce projection from %d - %d to %d - %d", this->subZ, this->subZ + pGeom->getDetectorRowCount(), this->subZ + _vmin, this->subZ + _vmax); + + return sub; +} + + +CCompositeGeometryManager::TPartList CCompositeGeometryManager::CProjectionPart::split(size_t maxSize, int div) +{ + TPartList ret; + + if (true) { + // Split in vertical direction only at first, until we figure out + // a model for splitting in other directions + + size_t sliceSize = ((size_t) pGeom->getDetectorColCount()) * pGeom->getProjectionCount(); + int sliceCount = pGeom->getDetectorRowCount(); + size_t blockSize = computeVerticalSplit(maxSize / sliceSize, div, sliceCount); + + int rem = sliceCount % blockSize; + + for (int z = -(rem / 2); z < sliceCount; z += blockSize) { + int newsubZ = z; + if (newsubZ < 0) newsubZ = 0; + int endZ = z + blockSize; + if (endZ > sliceCount) endZ = sliceCount; + int size = endZ - newsubZ; + + CProjectionPart *sub = new CProjectionPart(); + sub->subX = this->subX; + sub->subY = this->subY; + sub->subZ = this->subZ + newsubZ; + + ASTRA_DEBUG("ProjectionPart split %d %d %d -> %p", sub->subX, sub->subY, sub->subZ, (void*)sub); + + sub->pData = pData; + + sub->pGeom = getSubProjectionGeometry(pGeom, newsubZ, size); + + ret.push_back(boost::shared_ptr<CPart>(sub)); + } + } + + return ret; + +} + +CCompositeGeometryManager::CProjectionPart* CCompositeGeometryManager::CProjectionPart::clone() const +{ + return new CProjectionPart(*this); +} + + +bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, + CFloat32ProjectionData3DMemory *pProjData) +{ + ASTRA_DEBUG("CCompositeGeometryManager::doFP"); + // Create single job for FP + // Run result + + CVolumePart *input = new CVolumePart(); + input->pData = pVolData; + input->subX = 0; + input->subY = 0; + input->subZ = 0; + input->pGeom = pVolData->getGeometry()->clone(); + ASTRA_DEBUG("Main FP VolumePart -> %p", (void*)input); + + CProjectionPart *output = new CProjectionPart(); + output->pData = pProjData; + output->subX = 0; + output->subY = 0; + output->subZ = 0; + output->pGeom = pProjData->getGeometry()->clone(); + ASTRA_DEBUG("Main FP ProjectionPart -> %p", (void*)output); + + SJob FP; + FP.pInput = boost::shared_ptr<CPart>(input); + FP.pOutput = boost::shared_ptr<CPart>(output); + FP.pProjector = pProjector; + FP.eType = SJob::JOB_FP; + FP.eMode = SJob::MODE_SET; + + TJobList L; + L.push_back(FP); + + return doJobs(L); +} + +bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, + CFloat32ProjectionData3DMemory *pProjData) +{ + ASTRA_DEBUG("CCompositeGeometryManager::doBP"); + // Create single job for BP + // Run result + + CProjectionPart *input = new CProjectionPart(); + input->pData = pProjData; + input->subX = 0; + input->subY = 0; + input->subZ = 0; + input->pGeom = pProjData->getGeometry()->clone(); + + CVolumePart *output = new CVolumePart(); + output->pData = pVolData; + output->subX = 0; + output->subY = 0; + output->subZ = 0; + output->pGeom = pVolData->getGeometry()->clone(); + + SJob BP; + BP.pInput = boost::shared_ptr<CPart>(input); + BP.pOutput = boost::shared_ptr<CPart>(output); + BP.pProjector = pProjector; + BP.eType = SJob::JOB_BP; + BP.eMode = SJob::MODE_SET; + + TJobList L; + L.push_back(BP); + + return doJobs(L); +} + +bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData) +{ + ASTRA_DEBUG("CCompositeGeometryManager::doFP, multi-volume"); + + std::vector<CFloat32VolumeData3DMemory *>::const_iterator i; + std::vector<boost::shared_ptr<CPart> > inputs; + + for (i = volData.begin(); i != volData.end(); ++i) { + CVolumePart *input = new CVolumePart(); + input->pData = *i; + input->subX = 0; + input->subY = 0; + input->subZ = 0; + input->pGeom = (*i)->getGeometry()->clone(); + + inputs.push_back(boost::shared_ptr<CPart>(input)); + } + + std::vector<CFloat32ProjectionData3DMemory *>::const_iterator j; + std::vector<boost::shared_ptr<CPart> > outputs; + + for (j = projData.begin(); j != projData.end(); ++j) { + CProjectionPart *output = new CProjectionPart(); + output->pData = *j; + output->subX = 0; + output->subY = 0; + output->subZ = 0; + output->pGeom = (*j)->getGeometry()->clone(); + + outputs.push_back(boost::shared_ptr<CPart>(output)); + } + + std::vector<boost::shared_ptr<CPart> >::iterator i2; + std::vector<boost::shared_ptr<CPart> >::iterator j2; + TJobList L; + + for (i2 = outputs.begin(); i2 != outputs.end(); ++i2) { + SJob FP; + FP.eMode = SJob::MODE_SET; + for (j2 = inputs.begin(); j2 != inputs.end(); ++j2) { + FP.pInput = *j2; + FP.pOutput = *i2; + FP.pProjector = pProjector; + FP.eType = SJob::JOB_FP; + L.push_back(FP); + + // Set first, add rest + FP.eMode = SJob::MODE_ADD; + } + } + + return doJobs(L); +} + +bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData) +{ + ASTRA_DEBUG("CCompositeGeometryManager::doBP, multi-volume"); + + + std::vector<CFloat32VolumeData3DMemory *>::const_iterator i; + std::vector<boost::shared_ptr<CPart> > outputs; + + for (i = volData.begin(); i != volData.end(); ++i) { + CVolumePart *output = new CVolumePart(); + output->pData = *i; + output->subX = 0; + output->subY = 0; + output->subZ = 0; + output->pGeom = (*i)->getGeometry()->clone(); + + outputs.push_back(boost::shared_ptr<CPart>(output)); + } + + std::vector<CFloat32ProjectionData3DMemory *>::const_iterator j; + std::vector<boost::shared_ptr<CPart> > inputs; + + for (j = projData.begin(); j != projData.end(); ++j) { + CProjectionPart *input = new CProjectionPart(); + input->pData = *j; + input->subX = 0; + input->subY = 0; + input->subZ = 0; + input->pGeom = (*j)->getGeometry()->clone(); + + inputs.push_back(boost::shared_ptr<CPart>(input)); + } + + std::vector<boost::shared_ptr<CPart> >::iterator i2; + std::vector<boost::shared_ptr<CPart> >::iterator j2; + TJobList L; + + for (i2 = outputs.begin(); i2 != outputs.end(); ++i2) { + SJob BP; + BP.eMode = SJob::MODE_SET; + for (j2 = inputs.begin(); j2 != inputs.end(); ++j2) { + BP.pInput = *j2; + BP.pOutput = *i2; + BP.pProjector = pProjector; + BP.eType = SJob::JOB_BP; + L.push_back(BP); + + // Set first, add rest + BP.eMode = SJob::MODE_ADD; + } + } + + return doJobs(L); +} + + + + +bool CCompositeGeometryManager::doJobs(TJobList &jobs) +{ + ASTRA_DEBUG("CCompositeGeometryManager::doJobs"); + + // Sort job list into job set by output part + TJobSet jobset; + + for (TJobList::iterator i = jobs.begin(); i != jobs.end(); ++i) { + jobset[i->pOutput.get()].push_back(*i); + } + + size_t maxSize = astraCUDA3d::availableGPUMemory(); + if (maxSize == 0) { + ASTRA_WARN("Unable to get available GPU memory. Defaulting to 1GB."); + maxSize = 1024 * 1024 * 1024; + } else { + ASTRA_DEBUG("Detected %lu bytes of GPU memory", maxSize); + } + maxSize = (maxSize * 9) / 10; + + maxSize /= sizeof(float); + int div = 1; + + // TODO: Multi-GPU support + + // Split jobs to fit + TJobSet split; + splitJobs(jobset, maxSize, div, split); + jobset.clear(); + + // Run jobs + + for (TJobSet::iterator iter = split.begin(); iter != split.end(); ++iter) { + + CPart* output = iter->first; + TJobList& L = iter->second; + + assert(!L.empty()); + + bool zero = L.begin()->eMode == SJob::MODE_SET; + + size_t outx, outy, outz; + output->getDims(outx, outy, outz); + + if (L.begin()->eType == SJob::JOB_NOP) { + // just zero output? + if (zero) { + for (size_t z = 0; z < outz; ++z) { + for (size_t y = 0; y < outy; ++y) { + float* ptr = output->pData->getData(); + ptr += (z + output->subX) * (size_t)output->pData->getHeight() * (size_t)output->pData->getWidth(); + ptr += (y + output->subY) * (size_t)output->pData->getWidth(); + ptr += output->subX; + memset(ptr, 0, sizeof(float) * outx); + } + } + } + continue; + } + + + astraCUDA3d::SSubDimensions3D dstdims; + dstdims.nx = output->pData->getWidth(); + dstdims.pitch = dstdims.nx; + dstdims.ny = output->pData->getHeight(); + dstdims.nz = output->pData->getDepth(); + dstdims.subnx = outx; + dstdims.subny = outy; + dstdims.subnz = outz; + ASTRA_DEBUG("dstdims: %d,%d,%d in %d,%d,%d", dstdims.subnx, dstdims.subny, dstdims.subnz, dstdims.nx, dstdims.ny, dstdims.nz); + dstdims.subx = output->subX; + dstdims.suby = output->subY; + dstdims.subz = output->subZ; + float *dst = output->pData->getData(); + + astraCUDA3d::MemHandle3D outputMem = astraCUDA3d::allocateGPUMemory(outx, outy, outz, zero ? astraCUDA3d::INIT_ZERO : astraCUDA3d::INIT_NO); + bool ok = outputMem; + + for (TJobList::iterator i = L.begin(); i != L.end(); ++i) { + SJob &j = *i; + + assert(j.pInput); + + CCudaProjector3D *projector = dynamic_cast<CCudaProjector3D*>(j.pProjector); + Cuda3DProjectionKernel projKernel = ker3d_default; + int detectorSuperSampling = 1; + int voxelSuperSampling = 1; + if (projector) { + projKernel = projector->getProjectionKernel(); + detectorSuperSampling = projector->getDetectorSuperSampling(); + voxelSuperSampling = projector->getVoxelSuperSampling(); + } + + size_t inx, iny, inz; + j.pInput->getDims(inx, iny, inz); + astraCUDA3d::MemHandle3D inputMem = astraCUDA3d::allocateGPUMemory(inx, iny, inz, astraCUDA3d::INIT_NO); + + astraCUDA3d::SSubDimensions3D srcdims; + srcdims.nx = j.pInput->pData->getWidth(); + srcdims.pitch = srcdims.nx; + srcdims.ny = j.pInput->pData->getHeight(); + srcdims.nz = j.pInput->pData->getDepth(); + srcdims.subnx = inx; + srcdims.subny = iny; + srcdims.subnz = inz; + srcdims.subx = j.pInput->subX; + srcdims.suby = j.pInput->subY; + srcdims.subz = j.pInput->subZ; + const float *src = j.pInput->pData->getDataConst(); + + ok = astraCUDA3d::copyToGPUMemory(src, inputMem, srcdims); + if (!ok) ASTRA_ERROR("Error copying input data to GPU"); + + if (j.eType == SJob::JOB_FP) { + assert(dynamic_cast<CVolumePart*>(j.pInput.get())); + assert(dynamic_cast<CProjectionPart*>(j.pOutput.get())); + + ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing FP"); + + ok = astraCUDA3d::FP(((CProjectionPart*)j.pOutput.get())->pGeom, outputMem, ((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 == SJob::JOB_BP) { + assert(dynamic_cast<CVolumePart*>(j.pOutput.get())); + assert(dynamic_cast<CProjectionPart*>(j.pInput.get())); + + ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing BP"); + + ok = astraCUDA3d::BP(((CProjectionPart*)j.pInput.get())->pGeom, inputMem, ((CVolumePart*)j.pOutput.get())->pGeom, outputMem, voxelSuperSampling); + if (!ok) ASTRA_ERROR("Error performing sub-BP"); + ASTRA_DEBUG("CCompositeGeometryManager::doJobs: BP done"); + } else { + assert(false); + } + + ok = astraCUDA3d::freeGPUMemory(inputMem); + if (!ok) ASTRA_ERROR("Error freeing GPU memory"); + + } + + ok = astraCUDA3d::copyFromGPUMemory(dst, outputMem, dstdims); + if (!ok) ASTRA_ERROR("Error copying output data from GPU"); + + ok = astraCUDA3d::freeGPUMemory(outputMem); + if (!ok) ASTRA_ERROR("Error freeing GPU memory"); + } + + return true; +} + + + +} + +#endif diff --git a/src/ConeProjectionGeometry3D.cpp b/src/ConeProjectionGeometry3D.cpp index dd22eba..18f0f8a 100644 --- a/src/ConeProjectionGeometry3D.cpp +++ b/src/ConeProjectionGeometry3D.cpp @@ -29,6 +29,7 @@ $Id$ #include "astra/ConeProjectionGeometry3D.h" #include "astra/Logging.h" +#include "astra/GeometryUtil3D.h" #include <boost/lexical_cast.hpp> #include <cstring> @@ -230,14 +231,14 @@ CVector3D CConeProjectionGeometry3D::getProjectionDirection(int _iProjectionInde return ret; } -void CConeProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, - int iAngleIndex, - float32 &fU, float32 &fV) const +void CConeProjectionGeometry3D::projectPoint(double fX, double fY, double fZ, + int iAngleIndex, + double &fU, double &fV) const { ASTRA_ASSERT(iAngleIndex >= 0); ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); - float alpha = m_pfProjectionAngles[iAngleIndex]; + double alpha = m_pfProjectionAngles[iAngleIndex]; // Project point onto optical axis @@ -245,14 +246,14 @@ void CConeProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, // Vector source->origin is (-sin(alpha), cos(alpha)) // Distance from source, projected on optical axis - float fD = -sin(alpha) * fX + cos(alpha) * fY + m_fOriginSourceDistance; + double fD = -sin(alpha) * fX + cos(alpha) * fY + m_fOriginSourceDistance; // Scale fZ to detector plane fV = detectorOffsetYToRowIndexFloat( (fZ * (m_fOriginSourceDistance + m_fOriginDetectorDistance)) / fD ); // Orthogonal distance in XY-plane to optical axis - float fS = cos(alpha) * fX + sin(alpha) * fY; + double fS = cos(alpha) * fX + sin(alpha) * fY; // Scale fS to detector plane fU = detectorOffsetXToColIndexFloat( (fS * (m_fOriginSourceDistance + m_fOriginDetectorDistance)) / fD ); @@ -261,5 +262,84 @@ void CConeProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, } +void CConeProjectionGeometry3D::backprojectPointX(int iAngleIndex, double fU, double fV, + double fX, double &fY, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SConeProjection *projs = genConeProjections(1, m_iDetectorColCount, m_iDetectorRowCount, + m_fOriginSourceDistance, + m_fOriginDetectorDistance, + m_fDetectorSpacingX, m_fDetectorSpacingY, + &m_pfProjectionAngles[iAngleIndex]); + + SConeProjection &proj = projs[0]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - proj.fSrcX) / (px - proj.fSrcX); + + fY = proj.fSrcY + a * (py - proj.fSrcY); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); + + delete[] projs; +} + +void CConeProjectionGeometry3D::backprojectPointY(int iAngleIndex, double fU, double fV, + double fY, double &fX, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SConeProjection *projs = genConeProjections(1, m_iDetectorColCount, m_iDetectorRowCount, + m_fOriginSourceDistance, + m_fOriginDetectorDistance, + m_fDetectorSpacingX, m_fDetectorSpacingY, + &m_pfProjectionAngles[iAngleIndex]); + + SConeProjection &proj = projs[0]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - proj.fSrcY) / (py - proj.fSrcY); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); + + delete[] projs; +} + +void CConeProjectionGeometry3D::backprojectPointZ(int iAngleIndex, double fU, double fV, + double fZ, double &fX, double &fY) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SConeProjection *projs = genConeProjections(1, m_iDetectorColCount, m_iDetectorRowCount, + m_fOriginSourceDistance, + m_fOriginDetectorDistance, + m_fDetectorSpacingX, m_fDetectorSpacingY, + &m_pfProjectionAngles[iAngleIndex]); + + SConeProjection &proj = projs[0]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - proj.fSrcZ) / (pz - proj.fSrcZ); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fY = proj.fSrcY + a * (py - proj.fSrcY); + + delete[] projs; +} + + } // end namespace astra diff --git a/src/ConeVecProjectionGeometry3D.cpp b/src/ConeVecProjectionGeometry3D.cpp index 47ed630..86e3bd6 100644 --- a/src/ConeVecProjectionGeometry3D.cpp +++ b/src/ConeVecProjectionGeometry3D.cpp @@ -241,9 +241,9 @@ CVector3D CConeVecProjectionGeometry3D::getProjectionDirection(int _iProjectionI return CVector3D(p.fDetSX + (u+0.5)*p.fDetUX + (v+0.5)*p.fDetVX - p.fSrcX, p.fDetSY + (u+0.5)*p.fDetUY + (v+0.5)*p.fDetVY - p.fSrcY, p.fDetSZ + (u+0.5)*p.fDetUZ + (v+0.5)*p.fDetVZ - p.fSrcZ); } -void CConeVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, +void CConeVecProjectionGeometry3D::projectPoint(double fX, double fY, double fZ, int iAngleIndex, - float32 &fU, float32 &fV) const + double &fU, double &fV) const { ASTRA_ASSERT(iAngleIndex >= 0); ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); @@ -262,6 +262,60 @@ void CConeVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 } +void CConeVecProjectionGeometry3D::backprojectPointX(int iAngleIndex, double fU, double fV, + double fX, double &fY, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SConeProjection &proj = m_pProjectionAngles[iAngleIndex]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - proj.fSrcX) / (px - proj.fSrcX); + + fY = proj.fSrcY + a * (py - proj.fSrcY); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); +} + +void CConeVecProjectionGeometry3D::backprojectPointY(int iAngleIndex, double fU, double fV, + double fY, double &fX, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SConeProjection &proj = m_pProjectionAngles[iAngleIndex]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - proj.fSrcY) / (py - proj.fSrcY); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); +} + +void CConeVecProjectionGeometry3D::backprojectPointZ(int iAngleIndex, double fU, double fV, + double fZ, double &fX, double &fY) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SConeProjection &proj = m_pProjectionAngles[iAngleIndex]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - proj.fSrcZ) / (pz - proj.fSrcZ); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fY = proj.fSrcY + a * (py - proj.fSrcY); +} + //---------------------------------------------------------------------------------------- bool CConeVecProjectionGeometry3D::_check() diff --git a/src/CudaBackProjectionAlgorithm.cpp b/src/CudaBackProjectionAlgorithm.cpp index 365e058..a73f895 100644 --- a/src/CudaBackProjectionAlgorithm.cpp +++ b/src/CudaBackProjectionAlgorithm.cpp @@ -76,10 +76,9 @@ bool CCudaBackProjectionAlgorithm::initialize(const Config& _cfg) // Initialize - C++ bool CCudaBackProjectionAlgorithm::initialize(CProjector2D* _pProjector, CFloat32ProjectionData2D* _pSinogram, - CFloat32VolumeData2D* _pReconstruction, - int _iGPUindex, int _iPixelSuperSampling) + CFloat32VolumeData2D* _pReconstruction) { - m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, 1, _iPixelSuperSampling); + m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction); if (!m_bIsInitialized) return false; diff --git a/src/CudaBackProjectionAlgorithm3D.cpp b/src/CudaBackProjectionAlgorithm3D.cpp index e8e0433..ce8e111 100644 --- a/src/CudaBackProjectionAlgorithm3D.cpp +++ b/src/CudaBackProjectionAlgorithm3D.cpp @@ -37,6 +37,9 @@ $Id$ #include "astra/ParallelProjectionGeometry3D.h" #include "astra/ParallelVecProjectionGeometry3D.h" #include "astra/ConeVecProjectionGeometry3D.h" +#include "astra/CompositeGeometryManager.h" + +#include "astra/Logging.h" #include "../cuda/3d/astra3d.h" @@ -87,6 +90,24 @@ bool CCudaBackProjectionAlgorithm3D::_check() } //--------------------------------------------------------------------------------------- +void CCudaBackProjectionAlgorithm3D::initializeFromProjector() +{ + m_iVoxelSuperSampling = 1; + m_iGPUIndex = -1; + + CCudaProjector3D* pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector3D passed to BP3D_CUDA"); + } + } else { + m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } + +} + +//--------------------------------------------------------------------------------------- // Initialize - Config bool CCudaBackProjectionAlgorithm3D::initialize(const Config& _cfg) { @@ -103,32 +124,21 @@ bool CCudaBackProjectionAlgorithm3D::initialize(const Config& _cfg) return false; } - CCudaProjector3D* pCudaProjector = 0; - pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector); - if (!pCudaProjector) { - // TODO: Report - } + initializeFromProjector(); - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - CC.markOptionParsed("GPUindex"); - - - m_iVoxelSuperSampling = 1; - if (pCudaProjector) - m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + // Deprecated options m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); CC.markOptionParsed("VoxelSuperSampling"); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); - CFloat32ProjectionData3DMemory* pSinoMem = dynamic_cast<CFloat32ProjectionData3DMemory*>(m_pSinogram); - ASTRA_ASSERT(pSinoMem); - const CProjectionGeometry3D* projgeom = pSinoMem->getGeometry(); -const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(projgeom); - const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom); - if (parvec3dgeom || par3dgeom) { - // This option is only supported for Par3D currently - m_bSIRTWeighting = _cfg.self.getOptionBool("SIRTWeighting", false); - CC.markOptionParsed("SIRTWeighting"); - } + + + m_bSIRTWeighting = _cfg.self.getOptionBool("SIRTWeighting", false); + CC.markOptionParsed("SIRTWeighting"); // success m_bIsInitialized = _check(); @@ -151,6 +161,8 @@ bool CCudaBackProjectionAlgorithm3D::initialize(CProjector3D* _pProjector, m_pSinogram = _pSinogram; m_pReconstruction = _pReconstruction; + initializeFromProjector(); + // success m_bIsInitialized = _check(); return m_bIsInitialized; @@ -184,88 +196,24 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations) ASTRA_ASSERT(pReconMem); const CProjectionGeometry3D* projgeom = pSinoMem->getGeometry(); - const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(projgeom); - const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(projgeom); - const CConeVecProjectionGeometry3D* conevecgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(projgeom); - const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom); const CVolumeGeometry3D& volgeom = *pReconMem->getGeometry(); - if (conegeom) { - astraCudaConeBP(pReconMem->getData(), pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - conegeom->getProjectionCount(), - conegeom->getDetectorColCount(), - conegeom->getDetectorRowCount(), - conegeom->getOriginSourceDistance(), - conegeom->getOriginDetectorDistance(), - conegeom->getDetectorSpacingX(), - conegeom->getDetectorSpacingY(), - conegeom->getProjectionAngles(), - m_iGPUIndex, m_iVoxelSuperSampling); - } else if (par3dgeom) { - if (!m_bSIRTWeighting) { - astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - par3dgeom->getProjectionCount(), - par3dgeom->getDetectorColCount(), - par3dgeom->getDetectorRowCount(), - par3dgeom->getDetectorSpacingX(), - par3dgeom->getDetectorSpacingY(), - par3dgeom->getProjectionAngles(), - m_iGPUIndex, m_iVoxelSuperSampling); - } else { - astraCudaPar3DBP_SIRTWeighted(pReconMem->getData(), - pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - par3dgeom->getProjectionCount(), - par3dgeom->getDetectorColCount(), - par3dgeom->getDetectorRowCount(), - par3dgeom->getDetectorSpacingX(), - par3dgeom->getDetectorSpacingY(), - par3dgeom->getProjectionAngles(), - m_iGPUIndex, m_iVoxelSuperSampling); - } - } else if (parvec3dgeom) { - if (!m_bSIRTWeighting) { - astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - parvec3dgeom->getProjectionCount(), - parvec3dgeom->getDetectorColCount(), - parvec3dgeom->getDetectorRowCount(), - parvec3dgeom->getProjectionVectors(), - m_iGPUIndex, m_iVoxelSuperSampling); - } else { - astraCudaPar3DBP_SIRTWeighted(pReconMem->getData(), - pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - parvec3dgeom->getProjectionCount(), - parvec3dgeom->getDetectorColCount(), - parvec3dgeom->getDetectorRowCount(), - parvec3dgeom->getProjectionVectors(), - m_iGPUIndex, m_iVoxelSuperSampling); - } - } else if (conevecgeom) { - astraCudaConeBP(pReconMem->getData(), pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - conevecgeom->getProjectionCount(), - conevecgeom->getDetectorColCount(), - conevecgeom->getDetectorRowCount(), - conevecgeom->getProjectionVectors(), - m_iGPUIndex, m_iVoxelSuperSampling); + if (m_bSIRTWeighting) { + astraCudaBP_SIRTWeighted(pReconMem->getData(), + pSinoMem->getDataConst(), + &volgeom, projgeom, + m_iGPUIndex, m_iVoxelSuperSampling); } else { - ASTRA_ASSERT(false); + +#if 1 + CCompositeGeometryManager cgm; + + cgm.doBP(m_pProjector, pReconMem, pSinoMem); +#else + astraCudaBP(pReconMem->getData(), pSinoMem->getDataConst(), + &volgeom, projgeom, + m_iGPUIndex, m_iVoxelSuperSampling); +#endif } } diff --git a/src/CudaCglsAlgorithm.cpp b/src/CudaCglsAlgorithm.cpp index 0cedff6..9dd4f78 100644 --- a/src/CudaCglsAlgorithm.cpp +++ b/src/CudaCglsAlgorithm.cpp @@ -77,11 +77,9 @@ bool CCudaCglsAlgorithm::initialize(const Config& _cfg) // Initialize - C++ bool CCudaCglsAlgorithm::initialize(CProjector2D* _pProjector, CFloat32ProjectionData2D* _pSinogram, - CFloat32VolumeData2D* _pReconstruction, - int _iGPUindex, int _iDetectorSuperSampling, - int _iPixelSuperSampling) + CFloat32VolumeData2D* _pReconstruction) { - m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, _iDetectorSuperSampling, _iPixelSuperSampling); + m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction); if (!m_bIsInitialized) return false; diff --git a/src/CudaCglsAlgorithm3D.cpp b/src/CudaCglsAlgorithm3D.cpp index f527dc5..abc18d1 100644 --- a/src/CudaCglsAlgorithm3D.cpp +++ b/src/CudaCglsAlgorithm3D.cpp @@ -37,6 +37,8 @@ $Id$ #include "astra/ParallelVecProjectionGeometry3D.h" #include "astra/ConeVecProjectionGeometry3D.h" +#include "astra/Logging.h" + #include "../cuda/3d/astra3d.h" using namespace std; @@ -90,6 +92,26 @@ bool CCudaCglsAlgorithm3D::_check() } //--------------------------------------------------------------------------------------- +void CCudaCglsAlgorithm3D::initializeFromProjector() +{ + m_iVoxelSuperSampling = 1; + m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; + + CCudaProjector3D* pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector3D passed to CGLS3D_CUDA"); + } + } else { + m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } + +} + +//--------------------------------------------------------------------------------------- // Initialize - Config bool CCudaCglsAlgorithm3D::initialize(const Config& _cfg) { @@ -107,27 +129,20 @@ bool CCudaCglsAlgorithm3D::initialize(const Config& _cfg) return false; } - CCudaProjector3D* pCudaProjector = 0; - pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector); - if (!pCudaProjector) { - // TODO: Report - } + initializeFromProjector(); - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - CC.markOptionParsed("GPUindex"); - - m_iVoxelSuperSampling = 1; - m_iDetectorSuperSampling = 1; - if (pCudaProjector) { - // New interface - m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); - m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); - } // Deprecated options m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling); m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); CC.markOptionParsed("VoxelSuperSampling"); CC.markOptionParsed("DetectorSuperSampling"); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); + + m_pCgls = new AstraCGLS3d(); @@ -155,6 +170,8 @@ bool CCudaCglsAlgorithm3D::initialize(CProjector3D* _pProjector, m_pSinogram = _pSinogram; m_pReconstruction = _pReconstruction; + initializeFromProjector(); + m_pCgls = new AstraCGLS3d; m_bAstraCGLSInit = false; @@ -187,9 +204,6 @@ void CCudaCglsAlgorithm3D::run(int _iNrIterations) ASTRA_ASSERT(m_bIsInitialized); const CProjectionGeometry3D* projgeom = m_pSinogram->getGeometry(); - const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(projgeom); - const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom); - const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(projgeom); const CVolumeGeometry3D& volgeom = *m_pReconstruction->getGeometry(); bool ok = true; @@ -198,41 +212,7 @@ void CCudaCglsAlgorithm3D::run(int _iNrIterations) ok &= m_pCgls->setGPUIndex(m_iGPUIndex); - ok &= m_pCgls->setReconstructionGeometry(volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount()); -/* - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - float fOriginSourceDistance, - float fOriginDetectorDistance, - float fDetUSize, - float fDetVSize, - const float *pfAngles) -*/ - if (conegeom) { - ok &= m_pCgls->setConeGeometry(conegeom->getProjectionCount(), - conegeom->getDetectorColCount(), - conegeom->getDetectorRowCount(), - conegeom->getOriginSourceDistance(), - conegeom->getOriginDetectorDistance(), - conegeom->getDetectorSpacingX(), - conegeom->getDetectorSpacingY(), - conegeom->getProjectionAngles()); - } else if (parvec3dgeom) { - ok &= m_pCgls->setPar3DGeometry(parvec3dgeom->getProjectionCount(), - parvec3dgeom->getDetectorColCount(), - parvec3dgeom->getDetectorRowCount(), - parvec3dgeom->getProjectionVectors()); - } else if (conevec3dgeom) { - ok &= m_pCgls->setConeGeometry(conevec3dgeom->getProjectionCount(), - conevec3dgeom->getDetectorColCount(), - conevec3dgeom->getDetectorRowCount(), - conevec3dgeom->getProjectionVectors()); - } else { - ASTRA_ASSERT(false); - } + ok &= m_pCgls->setGeometry(&volgeom, projgeom); ok &= m_pCgls->enableSuperSampling(m_iVoxelSuperSampling, m_iDetectorSuperSampling); diff --git a/src/CudaEMAlgorithm.cpp b/src/CudaEMAlgorithm.cpp index 5c71f3d..d0afd80 100644 --- a/src/CudaEMAlgorithm.cpp +++ b/src/CudaEMAlgorithm.cpp @@ -76,11 +76,9 @@ bool CCudaEMAlgorithm::initialize(const Config& _cfg) // Initialize - C++ bool CCudaEMAlgorithm::initialize(CProjector2D* _pProjector, CFloat32ProjectionData2D* _pSinogram, - CFloat32VolumeData2D* _pReconstruction, - int _iGPUindex, int _iDetectorSuperSampling, - int _iPixelSuperSampling) + CFloat32VolumeData2D* _pReconstruction) { - m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, _iDetectorSuperSampling, _iPixelSuperSampling); + m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction); if (!m_bIsInitialized) return false; diff --git a/src/CudaFDKAlgorithm3D.cpp b/src/CudaFDKAlgorithm3D.cpp index 667d926..1316daa 100644 --- a/src/CudaFDKAlgorithm3D.cpp +++ b/src/CudaFDKAlgorithm3D.cpp @@ -35,6 +35,8 @@ $Id$ #include "astra/CudaProjector3D.h" #include "astra/ConeProjectionGeometry3D.h" +#include "astra/Logging.h" + #include "../cuda/3d/astra3d.h" using namespace std; @@ -85,6 +87,24 @@ bool CCudaFDKAlgorithm3D::_check() } //--------------------------------------------------------------------------------------- +void CCudaFDKAlgorithm3D::initializeFromProjector() +{ + m_iVoxelSuperSampling = 1; + m_iGPUIndex = -1; + + CCudaProjector3D* pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector3D passed to FDK_CUDA"); + } + } else { + m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } + +} + +//--------------------------------------------------------------------------------------- // Initialize - Config bool CCudaFDKAlgorithm3D::initialize(const Config& _cfg) { @@ -101,20 +121,18 @@ bool CCudaFDKAlgorithm3D::initialize(const Config& _cfg) return false; } - CCudaProjector3D* pCudaProjector = 0; - pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector); - if (!pCudaProjector) { - // TODO: Report - } - - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - CC.markOptionParsed("GPUindex"); + initializeFromProjector(); - m_iVoxelSuperSampling = 1; - if (pCudaProjector) - m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + // Deprecated options m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); CC.markOptionParsed("VoxelSuperSampling"); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); + + m_bShortScan = _cfg.self.getOptionBool("ShortScan", false); CC.markOptionParsed("ShortScan"); @@ -182,17 +200,7 @@ void CCudaFDKAlgorithm3D::run(int _iNrIterations) bool ok = true; ok = astraCudaFDK(pReconMem->getData(), pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - conegeom->getProjectionCount(), - conegeom->getDetectorColCount(), - conegeom->getDetectorRowCount(), - conegeom->getOriginSourceDistance(), - conegeom->getOriginDetectorDistance(), - conegeom->getDetectorSpacingX(), - conegeom->getDetectorSpacingY(), - conegeom->getProjectionAngles(), + &volgeom, conegeom, m_bShortScan, m_iGPUIndex, m_iVoxelSuperSampling); ASTRA_ASSERT(ok); diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp index aac96d6..bcd70c4 100644 --- a/src/CudaFilteredBackProjectionAlgorithm.cpp +++ b/src/CudaFilteredBackProjectionAlgorithm.cpp @@ -67,6 +67,24 @@ CCudaFilteredBackProjectionAlgorithm::~CCudaFilteredBackProjectionAlgorithm() } } +void CCudaFilteredBackProjectionAlgorithm::initializeFromProjector() +{ + m_iPixelSuperSampling = 1; + m_iGPUIndex = -1; + + // Projector + CCudaProjector2D* pCudaProjector = dynamic_cast<CCudaProjector2D*>(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector2D passed to FBP_CUDA"); + } + } else { + m_iPixelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } + +} + bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) { ASTRA_ASSERT(_cfg.self); @@ -163,35 +181,30 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) } CC.markNodeParsed("FilterD"); // TODO: Only for some types! - // GPU number - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - CC.markOptionParsed("GPUindex"); - - m_iPixelSuperSampling = 1; - if (pCudaProjector) { - // New interface - m_iPixelSuperSampling = pCudaProjector->getVoxelSuperSampling(); - } - // Deprecated options - m_iPixelSuperSampling = (int)_cfg.self.getOptionNumerical("PixelSuperSampling", m_iPixelSuperSampling); - CC.markOptionParsed("PixelSuperSampling"); - - // Fan beam short scan mode if (m_pSinogram && dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pSinogram->getGeometry())) { m_bShortScan = (int)_cfg.self.getOptionBool("ShortScan", false); CC.markOptionParsed("ShortScan"); } + initializeFromProjector(); + + // Deprecated options + m_iPixelSuperSampling = (int)_cfg.self.getOptionNumerical("PixelSuperSampling", m_iPixelSuperSampling); + CC.markOptionParsed("PixelSuperSampling"); + // GPU number + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); m_pFBP = new AstraFBP; m_bAstraFBPInit = false; - // success - m_bIsInitialized = true; - return m_bIsInitialized; + return check(); } bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D * _pSinogram, CFloat32VolumeData2D * _pReconstruction, E_FBPFILTER _eFilter, const float * _pfFilter /* = NULL */, int _iFilterWidth /* = 0 */, int _iGPUIndex /* = 0 */, float _fFilterParameter /* = -1.0f */) @@ -241,7 +254,7 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D * m_fFilterParameter = _fFilterParameter; - return m_bIsInitialized; + return check(); } void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */) @@ -361,7 +374,7 @@ bool CCudaFilteredBackProjectionAlgorithm::check() ASTRA_CONFIG_CHECK(m_pReconstruction->isInitialized(), "FBP_CUDA", "Reconstruction Data Object Not Initialized."); // check gpu index - ASTRA_CONFIG_CHECK(m_iGPUIndex >= -1, "FBP_CUDA", "GPUIndex must be a non-negative integer."); + ASTRA_CONFIG_CHECK(m_iGPUIndex >= -1, "FBP_CUDA", "GPUIndex must be a non-negative integer or -1."); // check pixel supersampling ASTRA_CONFIG_CHECK(m_iPixelSuperSampling >= 0, "FBP_CUDA", "PixelSuperSampling must be a non-negative integer."); diff --git a/src/CudaForwardProjectionAlgorithm.cpp b/src/CudaForwardProjectionAlgorithm.cpp index b382f2e..9ca13ae 100644 --- a/src/CudaForwardProjectionAlgorithm.cpp +++ b/src/CudaForwardProjectionAlgorithm.cpp @@ -38,8 +38,11 @@ $Id$ #include <boost/lexical_cast.hpp> #include "astra/AstraObjectManager.h" +#include "astra/ParallelProjectionGeometry2D.h" #include "astra/FanFlatProjectionGeometry2D.h" #include "astra/FanFlatVecProjectionGeometry2D.h" +#include "astra/Float32ProjectionData2D.h" +#include "astra/Float32VolumeData2D.h" #include "astra/CudaProjector2D.h" #include "astra/Logging.h" @@ -66,6 +69,24 @@ CCudaForwardProjectionAlgorithm::~CCudaForwardProjectionAlgorithm() } //--------------------------------------------------------------------------------------- +void CCudaForwardProjectionAlgorithm::initializeFromProjector() +{ + m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; + + // Projector + CCudaProjector2D* pCudaProjector = dynamic_cast<CCudaProjector2D*>(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector2D passed to FP_CUDA"); + } + } else { + m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } +} + +//--------------------------------------------------------------------------------------- // Initialize - Config bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg) { @@ -74,14 +95,9 @@ bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg) // Projector XMLNode node = _cfg.self.getSingleNode("ProjectorId"); - CCudaProjector2D* pCudaProjector = 0; if (node) { int id = boost::lexical_cast<int>(node.getContent()); - CProjector2D *projector = CProjector2DManager::getSingleton().get(id); - pCudaProjector = dynamic_cast<CCudaProjector2D*>(projector); - if (!pCudaProjector) { - ASTRA_WARN("non-CUDA Projector2D passed to FP_CUDA"); - } + m_pProjector = CProjector2DManager::getSingleton().get(id); } CC.markNodeParsed("ProjectorId"); @@ -101,22 +117,18 @@ bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg) m_pVolume = dynamic_cast<CFloat32VolumeData2D*>(CData2DManager::getSingleton().get(id)); CC.markNodeParsed("VolumeDataId"); + initializeFromProjector(); + + // Deprecated options + m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling); + CC.markOptionParsed("DetectorSuperSampling"); // GPU number m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); - CC.markOptionParsed("GPUindex"); - if (!_cfg.self.hasOption("GPUindex")) - CC.markOptionParsed("GPUIndex"); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); - // Detector supersampling factor - m_iDetectorSuperSampling = 1; - if (pCudaProjector) { - // New interface - m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); - } - // Deprecated option - m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling); - CC.markOptionParsed("DetectorSuperSampling"); // return success @@ -125,20 +137,16 @@ bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg) //---------------------------------------------------------------------------------------- // Initialize - C++ -bool CCudaForwardProjectionAlgorithm::initialize(CProjectionGeometry2D* _pProjectionGeometry, - CVolumeGeometry2D* _pReconstructionGeometry, +bool CCudaForwardProjectionAlgorithm::initialize(CProjector2D* _pProjector, CFloat32VolumeData2D* _pVolume, - CFloat32ProjectionData2D* _pSinogram, - int _iGPUindex, int _iDetectorSuperSampling) + CFloat32ProjectionData2D* _pSinogram) { // store classes - //m_pProjectionGeometry = _pProjectionGeometry; - //m_pReconstructionGeometry = _pReconstructionGeometry; + m_pProjector = _pProjector; m_pVolume = _pVolume; m_pSinogram = _pSinogram; - m_iDetectorSuperSampling = _iDetectorSuperSampling; - m_iGPUIndex = _iGPUindex; + initializeFromProjector(); // return success return check(); diff --git a/src/CudaForwardProjectionAlgorithm3D.cpp b/src/CudaForwardProjectionAlgorithm3D.cpp index 46dab12..209f5a5 100644 --- a/src/CudaForwardProjectionAlgorithm3D.cpp +++ b/src/CudaForwardProjectionAlgorithm3D.cpp @@ -40,6 +40,8 @@ $Id$ #include "astra/ParallelVecProjectionGeometry3D.h" #include "astra/ConeVecProjectionGeometry3D.h" +#include "astra/CompositeGeometryManager.h" + #include "astra/Logging.h" #include "../cuda/3d/astra3d.h" @@ -72,6 +74,23 @@ CCudaForwardProjectionAlgorithm3D::~CCudaForwardProjectionAlgorithm3D() } //--------------------------------------------------------------------------------------- +void CCudaForwardProjectionAlgorithm3D::initializeFromProjector() +{ + m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; + + CCudaProjector3D* pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector3D passed to FP3D_CUDA"); + } + } else { + m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } +} + +//--------------------------------------------------------------------------------------- // Initialize - Config bool CCudaForwardProjectionAlgorithm3D::initialize(const Config& _cfg) { @@ -97,29 +116,21 @@ bool CCudaForwardProjectionAlgorithm3D::initialize(const Config& _cfg) // optional: projector node = _cfg.self.getSingleNode("ProjectorId"); - CCudaProjector3D* pCudaProjector = 0; m_pProjector = 0; if (node) { id = boost::lexical_cast<int>(node.getContent()); m_pProjector = CProjector3DManager::getSingleton().get(id); - pCudaProjector = dynamic_cast<CCudaProjector3D*>(CProjector3DManager::getSingleton().get(id)); - m_pProjector = pCudaProjector; - if (!pCudaProjector) { - // TODO: Report - } } CC.markNodeParsed("ProjectorId"); - // GPU number - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - CC.markOptionParsed("GPUindex"); - + initializeFromProjector(); - m_iDetectorSuperSampling = 1; - if (pCudaProjector) - m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); + // Deprecated options m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex); CC.markOptionParsed("DetectorSuperSampling"); + CC.markOptionParsed("GPUindex"); + // success m_bIsInitialized = check(); @@ -142,8 +153,15 @@ bool CCudaForwardProjectionAlgorithm3D::initialize(CProjector3D* _pProjector, m_pProjections = _pProjections; m_pVolume = _pVolume; - m_iDetectorSuperSampling = _iDetectorSuperSampling; - m_iGPUIndex = _iGPUindex; + CCudaProjector3D* pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector); + if (!pCudaProjector) { + // TODO: Report + m_iDetectorSuperSampling = _iDetectorSuperSampling; + m_iGPUIndex = _iGPUindex; + } else { + m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } // success m_bIsInitialized = check(); @@ -247,11 +265,13 @@ void CCudaForwardProjectionAlgorithm3D::run(int) // check initialized assert(m_bIsInitialized); +#if 1 + CCompositeGeometryManager cgm; + + cgm.doFP(m_pProjector, m_pVolume, m_pProjections); + +#else const CProjectionGeometry3D* projgeom = m_pProjections->getGeometry(); - const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(projgeom); - const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(projgeom); - const CConeVecProjectionGeometry3D* conevecgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(projgeom); - const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom); const CVolumeGeometry3D& volgeom = *m_pVolume->getGeometry(); Cuda3DProjectionKernel projKernel = ker3d_default; @@ -279,58 +299,10 @@ void CCudaForwardProjectionAlgorithm3D::run(int) } #endif - if (conegeom) { - astraCudaConeFP(m_pVolume->getDataConst(), m_pProjections->getData(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - conegeom->getProjectionCount(), - conegeom->getDetectorColCount(), - conegeom->getDetectorRowCount(), - conegeom->getOriginSourceDistance(), - conegeom->getOriginDetectorDistance(), - conegeom->getDetectorSpacingX(), - conegeom->getDetectorSpacingY(), - conegeom->getProjectionAngles(), - m_iGPUIndex, m_iDetectorSuperSampling); - } else if (par3dgeom) { - astraCudaPar3DFP(m_pVolume->getDataConst(), m_pProjections->getData(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - par3dgeom->getProjectionCount(), - par3dgeom->getDetectorColCount(), - par3dgeom->getDetectorRowCount(), - par3dgeom->getDetectorSpacingX(), - par3dgeom->getDetectorSpacingY(), - par3dgeom->getProjectionAngles(), - m_iGPUIndex, m_iDetectorSuperSampling, - projKernel); - } else if (parvec3dgeom) { - astraCudaPar3DFP(m_pVolume->getDataConst(), m_pProjections->getData(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - parvec3dgeom->getProjectionCount(), - parvec3dgeom->getDetectorColCount(), - parvec3dgeom->getDetectorRowCount(), - parvec3dgeom->getProjectionVectors(), - m_iGPUIndex, m_iDetectorSuperSampling, - projKernel); - } else if (conevecgeom) { - astraCudaConeFP(m_pVolume->getDataConst(), m_pProjections->getData(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - conevecgeom->getProjectionCount(), - conevecgeom->getDetectorColCount(), - conevecgeom->getDetectorRowCount(), - conevecgeom->getProjectionVectors(), - m_iGPUIndex, m_iDetectorSuperSampling); - } else { - ASTRA_ASSERT(false); - } - + astraCudaFP(m_pVolume->getDataConst(), m_pProjections->getData(), + &volgeom, projgeom, + m_iGPUIndex, m_iDetectorSuperSampling, projKernel); +#endif } diff --git a/src/CudaProjector2D.cpp b/src/CudaProjector2D.cpp index a26e32d..acf6000 100644 --- a/src/CudaProjector2D.cpp +++ b/src/CudaProjector2D.cpp @@ -61,6 +61,7 @@ void CCudaProjector2D::_clear() m_projectionKernel = ker2d_default; m_iVoxelSuperSampling = 1; m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; } //---------------------------------------------------------------------------------------- @@ -125,18 +126,18 @@ bool CCudaProjector2D::initialize(const Config& _cfg) m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", 1); CC.markOptionParsed("DetectorSuperSampling"); + // GPU number + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); + + m_bIsInitialized = _check(); return m_bIsInitialized; } -/* -bool CProjector2D::initialize(astra::CProjectionGeometry2D *, astra::CVolumeGeometry2D *) -{ - ASTRA_ASSERT(false); - - return false; -} -*/ std::string CCudaProjector2D::description() const { diff --git a/src/CudaProjector3D.cpp b/src/CudaProjector3D.cpp index d2fd74c..bbfbd34 100644 --- a/src/CudaProjector3D.cpp +++ b/src/CudaProjector3D.cpp @@ -64,6 +64,7 @@ void CCudaProjector3D::_clear() m_projectionKernel = ker3d_default; m_iVoxelSuperSampling = 1; m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; } //---------------------------------------------------------------------------------------- @@ -128,6 +129,12 @@ bool CCudaProjector3D::initialize(const Config& _cfg) m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", 1); CC.markOptionParsed("DetectorSuperSampling"); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); + m_bIsInitialized = _check(); return m_bIsInitialized; } diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp index 71b6637..71dddf7 100644 --- a/src/CudaReconstructionAlgorithm2D.cpp +++ b/src/CudaReconstructionAlgorithm2D.cpp @@ -84,111 +84,51 @@ void CCudaReconstructionAlgorithm2D::_clear() } //--------------------------------------------------------------------------------------- -// Initialize - Config -bool CCudaReconstructionAlgorithm2D::initialize(const Config& _cfg) +void CCudaReconstructionAlgorithm2D::initializeFromProjector() { - ASTRA_ASSERT(_cfg.self); - ConfigStackCheck<CAlgorithm> CC("CudaReconstructionAlgorithm2D", this, _cfg); - - // if already initialized, clear first - if (m_bIsInitialized) { - clear(); - } + m_iPixelSuperSampling = 1; + m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; // Projector - XMLNode node = _cfg.self.getSingleNode("ProjectorId"); - CCudaProjector2D* pCudaProjector = 0; - if (node) { - int id = boost::lexical_cast<int>(node.getContent()); - CProjector2D *projector = CProjector2DManager::getSingleton().get(id); - pCudaProjector = dynamic_cast<CCudaProjector2D*>(projector); - if (!pCudaProjector) { + CCudaProjector2D* pCudaProjector = dynamic_cast<CCudaProjector2D*>(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { ASTRA_WARN("non-CUDA Projector2D passed"); } - } - CC.markNodeParsed("ProjectorId"); - - - // sinogram data - node = _cfg.self.getSingleNode("ProjectionDataId"); - ASTRA_CONFIG_CHECK(node, "CudaSirt2", "No ProjectionDataId tag specified."); - int id = boost::lexical_cast<int>(node.getContent()); - m_pSinogram = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id)); - CC.markNodeParsed("ProjectionDataId"); - - // reconstruction data - node = _cfg.self.getSingleNode("ReconstructionDataId"); - ASTRA_CONFIG_CHECK(node, "CudaSirt2", "No ReconstructionDataId tag specified."); - id = boost::lexical_cast<int>(node.getContent()); - m_pReconstruction = dynamic_cast<CFloat32VolumeData2D*>(CData2DManager::getSingleton().get(id)); - CC.markNodeParsed("ReconstructionDataId"); - - // fixed mask - if (_cfg.self.hasOption("ReconstructionMaskId")) { - m_bUseReconstructionMask = true; - id = boost::lexical_cast<int>(_cfg.self.getOption("ReconstructionMaskId")); - m_pReconstructionMask = dynamic_cast<CFloat32VolumeData2D*>(CData2DManager::getSingleton().get(id)); - ASTRA_CONFIG_CHECK(m_pReconstructionMask, "CudaReconstruction2D", "Invalid ReconstructionMaskId."); - } - CC.markOptionParsed("ReconstructionMaskId"); - // fixed mask - if (_cfg.self.hasOption("SinogramMaskId")) { - m_bUseSinogramMask = true; - id = boost::lexical_cast<int>(_cfg.self.getOption("SinogramMaskId")); - m_pSinogramMask = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id)); - ASTRA_CONFIG_CHECK(m_pSinogramMask, "CudaReconstruction2D", "Invalid SinogramMaskId."); - } - CC.markOptionParsed("SinogramMaskId"); - - // Constraints - NEW - if (_cfg.self.hasOption("MinConstraint")) { - m_bUseMinConstraint = true; - m_fMinValue = _cfg.self.getOptionNumerical("MinConstraint", 0.0f); - CC.markOptionParsed("MinConstraint"); - } else { - // Constraint - OLD - m_bUseMinConstraint = _cfg.self.getOptionBool("UseMinConstraint", false); - CC.markOptionParsed("UseMinConstraint"); - if (m_bUseMinConstraint) { - m_fMinValue = _cfg.self.getOptionNumerical("MinConstraintValue", 0.0f); - CC.markOptionParsed("MinConstraintValue"); - } - } - if (_cfg.self.hasOption("MaxConstraint")) { - m_bUseMaxConstraint = true; - m_fMaxValue = _cfg.self.getOptionNumerical("MaxConstraint", 255.0f); - CC.markOptionParsed("MaxConstraint"); } else { - // Constraint - OLD - m_bUseMaxConstraint = _cfg.self.getOptionBool("UseMaxConstraint", false); - CC.markOptionParsed("UseMaxConstraint"); - if (m_bUseMaxConstraint) { - m_fMaxValue = _cfg.self.getOptionNumerical("MaxConstraintValue", 0.0f); - CC.markOptionParsed("MaxConstraintValue"); - } - } - - // GPU number - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); - CC.markOptionParsed("GPUindex"); - if (!_cfg.self.hasOption("GPUindex")) - CC.markOptionParsed("GPUIndex"); - - // Supersampling factors - m_iDetectorSuperSampling = 1; - m_iPixelSuperSampling = 1; - if (pCudaProjector) { - // New interface m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); m_iPixelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); } +} + +//--------------------------------------------------------------------------------------- +// Initialize - Config +bool CCudaReconstructionAlgorithm2D::initialize(const Config& _cfg) +{ + ASTRA_ASSERT(_cfg.self); + ConfigStackCheck<CAlgorithm> CC("CudaReconstructionAlgorithm2D", this, _cfg); + + m_bIsInitialized = CReconstructionAlgorithm2D::initialize(_cfg); + + if (!m_bIsInitialized) + return false; + + initializeFromProjector(); + // Deprecated options m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling); m_iPixelSuperSampling = (int)_cfg.self.getOptionNumerical("PixelSuperSampling", m_iPixelSuperSampling); CC.markOptionParsed("DetectorSuperSampling"); CC.markOptionParsed("PixelSuperSampling"); + // GPU number + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); return _check(); } @@ -199,32 +139,18 @@ bool CCudaReconstructionAlgorithm2D::initialize(CProjector2D* _pProjector, CFloat32ProjectionData2D* _pSinogram, CFloat32VolumeData2D* _pReconstruction) { - return initialize(_pProjector, _pSinogram, _pReconstruction, 0, 1); -} - -//--------------------------------------------------------------------------------------- -// Initialize - C++ -bool CCudaReconstructionAlgorithm2D::initialize(CProjector2D* _pProjector, - CFloat32ProjectionData2D* _pSinogram, - CFloat32VolumeData2D* _pReconstruction, - int _iGPUindex, - int _iDetectorSuperSampling, - int _iPixelSuperSampling) -{ // if already initialized, clear first if (m_bIsInitialized) { clear(); } - m_pProjector = 0; + m_pProjector = _pProjector; // required classes m_pSinogram = _pSinogram; m_pReconstruction = _pReconstruction; - m_iDetectorSuperSampling = _iDetectorSuperSampling; - m_iPixelSuperSampling = _iPixelSuperSampling; - m_iGPUIndex = _iGPUindex; + initializeFromProjector(); return _check(); } @@ -234,40 +160,13 @@ bool CCudaReconstructionAlgorithm2D::initialize(CProjector2D* _pProjector, // Check bool CCudaReconstructionAlgorithm2D::_check() { - // TODO: CLEAN UP - - - // check pointers - //ASTRA_CONFIG_CHECK(m_pProjector, "Reconstruction2D", "Invalid Projector Object."); - ASTRA_CONFIG_CHECK(m_pSinogram, "SIRT_CUDA", "Invalid Projection Data Object."); - ASTRA_CONFIG_CHECK(m_pReconstruction, "SIRT_CUDA", "Invalid Reconstruction Data Object."); - - // check initializations - //ASTRA_CONFIG_CHECK(m_pProjector->isInitialized(), "Reconstruction2D", "Projector Object Not Initialized."); - ASTRA_CONFIG_CHECK(m_pSinogram->isInitialized(), "SIRT_CUDA", "Projection Data Object Not Initialized."); - ASTRA_CONFIG_CHECK(m_pReconstruction->isInitialized(), "SIRT_CUDA", "Reconstruction Data Object Not Initialized."); - - ASTRA_CONFIG_CHECK(m_iDetectorSuperSampling >= 1, "SIRT_CUDA", "DetectorSuperSampling must be a positive integer."); - ASTRA_CONFIG_CHECK(m_iPixelSuperSampling >= 1, "SIRT_CUDA", "PixelSuperSampling must be a positive integer."); - ASTRA_CONFIG_CHECK(m_iGPUIndex >= -1, "SIRT_CUDA", "GPUIndex must be a non-negative integer."); - - // check compatibility between projector and data classes -// ASTRA_CONFIG_CHECK(m_pSinogram->getGeometry()->isEqual(m_pProjector->getProjectionGeometry()), "SIRT_CUDA", "Projection Data not compatible with the specified Projector."); -// ASTRA_CONFIG_CHECK(m_pReconstruction->getGeometry()->isEqual(m_pProjector->getVolumeGeometry()), "SIRT_CUDA", "Reconstruction Data not compatible with the specified Projector."); - - // todo: turn some of these back on + if (!CReconstructionAlgorithm2D::_check()) + return false; -// ASTRA_CONFIG_CHECK(m_pProjectionGeometry, "SIRT_CUDA", "ProjectionGeometry not specified."); -// ASTRA_CONFIG_CHECK(m_pProjectionGeometry->isInitialized(), "SIRT_CUDA", "ProjectionGeometry not initialized."); -// ASTRA_CONFIG_CHECK(m_pReconstructionGeometry, "SIRT_CUDA", "ReconstructionGeometry not specified."); -// ASTRA_CONFIG_CHECK(m_pReconstructionGeometry->isInitialized(), "SIRT_CUDA", "ReconstructionGeometry not initialized."); + ASTRA_CONFIG_CHECK(m_iDetectorSuperSampling >= 1, "CudaReconstructionAlgorithm2D", "DetectorSuperSampling must be a positive integer."); + ASTRA_CONFIG_CHECK(m_iPixelSuperSampling >= 1, "CudaReconstructionAlgorithm2D", "PixelSuperSampling must be a positive integer."); + ASTRA_CONFIG_CHECK(m_iGPUIndex >= -1, "CudaReconstructionAlgorithm2D", "GPUIndex must be a non-negative integer or -1."); - // check dimensions - //ASTRA_CONFIG_CHECK(m_pSinogram->getAngleCount() == m_pProjectionGeometry->getProjectionAngleCount(), "SIRT_CUDA", "Sinogram data object size mismatch."); - //ASTRA_CONFIG_CHECK(m_pSinogram->getDetectorCount() == m_pProjectionGeometry->getDetectorCount(), "SIRT_CUDA", "Sinogram data object size mismatch."); - //ASTRA_CONFIG_CHECK(m_pReconstruction->getWidth() == m_pReconstructionGeometry->getGridColCount(), "SIRT_CUDA", "Reconstruction data object size mismatch."); - //ASTRA_CONFIG_CHECK(m_pReconstruction->getHeight() == m_pReconstructionGeometry->getGridRowCount(), "SIRT_CUDA", "Reconstruction data object size mismatch."); - // check restrictions // TODO: check restrictions built into cuda code @@ -462,10 +361,18 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations) ASTRA_ASSERT(ok); - if (m_bUseMinConstraint) - ok &= m_pAlgo->setMinConstraint(m_fMinValue); - if (m_bUseMaxConstraint) - ok &= m_pAlgo->setMaxConstraint(m_fMaxValue); + if (m_bUseMinConstraint) { + bool ret = m_pAlgo->setMinConstraint(m_fMinValue); + if (!ret) { + ASTRA_WARN("This algorithm ignores MinConstraint"); + } + } + if (m_bUseMaxConstraint) { + bool ret= m_pAlgo->setMaxConstraint(m_fMaxValue); + if (!ret) { + ASTRA_WARN("This algorithm ignores MaxConstraint"); + } + } ok &= m_pAlgo->iterate(_iNrIterations); ASTRA_ASSERT(ok); diff --git a/src/CudaSartAlgorithm.cpp b/src/CudaSartAlgorithm.cpp index 8c0c6d7..d202847 100644 --- a/src/CudaSartAlgorithm.cpp +++ b/src/CudaSartAlgorithm.cpp @@ -116,10 +116,9 @@ bool CCudaSartAlgorithm::initialize(const Config& _cfg) // Initialize - C++ bool CCudaSartAlgorithm::initialize(CProjector2D* _pProjector, CFloat32ProjectionData2D* _pSinogram, - CFloat32VolumeData2D* _pReconstruction, - int _iGPUindex, int _iDetectorSuperSampling) + CFloat32VolumeData2D* _pReconstruction) { - m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, _iDetectorSuperSampling, 1); + m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction); if (!m_bIsInitialized) return false; diff --git a/src/CudaSirtAlgorithm.cpp b/src/CudaSirtAlgorithm.cpp index d424915..ab0a418 100644 --- a/src/CudaSirtAlgorithm.cpp +++ b/src/CudaSirtAlgorithm.cpp @@ -98,11 +98,9 @@ bool CCudaSirtAlgorithm::initialize(const Config& _cfg) // Initialize - C++ bool CCudaSirtAlgorithm::initialize(CProjector2D* _pProjector, CFloat32ProjectionData2D* _pSinogram, - CFloat32VolumeData2D* _pReconstruction, - int _iGPUindex, int _iDetectorSuperSampling, - int _iPixelSuperSampling) + CFloat32VolumeData2D* _pReconstruction) { - m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, _iDetectorSuperSampling, _iPixelSuperSampling); + m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction); if (!m_bIsInitialized) return false; diff --git a/src/CudaSirtAlgorithm3D.cpp b/src/CudaSirtAlgorithm3D.cpp index abbb9fd..1fa0da2 100644 --- a/src/CudaSirtAlgorithm3D.cpp +++ b/src/CudaSirtAlgorithm3D.cpp @@ -38,6 +38,8 @@ $Id$ #include "astra/ConeVecProjectionGeometry3D.h" #include "astra/CudaProjector3D.h" +#include "astra/Logging.h" + #include "../cuda/3d/astra3d.h" using namespace std; @@ -90,7 +92,27 @@ bool CCudaSirtAlgorithm3D::_check() return true; } -//--------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------------- +void CCudaSirtAlgorithm3D::initializeFromProjector() +{ + m_iVoxelSuperSampling = 1; + m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; + + CCudaProjector3D* pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector3D passed to SIRT3D_CUDA"); + } + } else { + m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } + +} + +//-------------------------------------------------------------------------------------- // Initialize - Config bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg) { @@ -108,28 +130,20 @@ bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg) return false; } - CCudaProjector3D* pCudaProjector = 0; - pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector); - if (!pCudaProjector) { - // TODO: Report - } + initializeFromProjector(); - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - CC.markOptionParsed("GPUindex"); - - - m_iVoxelSuperSampling = 1; - m_iDetectorSuperSampling = 1; - if (pCudaProjector) { - // New interface - m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); - m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); - } // Deprecated options m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling); m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); CC.markOptionParsed("VoxelSuperSampling"); CC.markOptionParsed("DetectorSuperSampling"); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); + + m_pSirt = new AstraSIRT3d(); @@ -189,10 +203,6 @@ void CCudaSirtAlgorithm3D::run(int _iNrIterations) ASTRA_ASSERT(m_bIsInitialized); const CProjectionGeometry3D* projgeom = m_pSinogram->getGeometry(); - const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(projgeom); - const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(projgeom); - const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom); - const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(projgeom); const CVolumeGeometry3D& volgeom = *m_pReconstruction->getGeometry(); bool ok = true; @@ -201,39 +211,7 @@ void CCudaSirtAlgorithm3D::run(int _iNrIterations) ok &= m_pSirt->setGPUIndex(m_iGPUIndex); - ok &= m_pSirt->setReconstructionGeometry(volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount()); - - if (conegeom) { - ok &= m_pSirt->setConeGeometry(conegeom->getProjectionCount(), - conegeom->getDetectorColCount(), - conegeom->getDetectorRowCount(), - conegeom->getOriginSourceDistance(), - conegeom->getOriginDetectorDistance(), - conegeom->getDetectorSpacingX(), - conegeom->getDetectorSpacingY(), - conegeom->getProjectionAngles()); - } else if (par3dgeom) { - ok &= m_pSirt->setPar3DGeometry(par3dgeom->getProjectionCount(), - par3dgeom->getDetectorColCount(), - par3dgeom->getDetectorRowCount(), - par3dgeom->getDetectorSpacingX(), - par3dgeom->getDetectorSpacingY(), - par3dgeom->getProjectionAngles()); - } else if (parvec3dgeom) { - ok &= m_pSirt->setPar3DGeometry(parvec3dgeom->getProjectionCount(), - parvec3dgeom->getDetectorColCount(), - parvec3dgeom->getDetectorRowCount(), - parvec3dgeom->getProjectionVectors()); - } else if (conevec3dgeom) { - ok &= m_pSirt->setConeGeometry(conevec3dgeom->getProjectionCount(), - conevec3dgeom->getDetectorColCount(), - conevec3dgeom->getDetectorRowCount(), - conevec3dgeom->getProjectionVectors()); - } else { - ASTRA_ASSERT(false); - } + ok &= m_pSirt->setGeometry(&volgeom, projgeom); ok &= m_pSirt->enableSuperSampling(m_iVoxelSuperSampling, m_iDetectorSuperSampling); diff --git a/src/Float32ProjectionData3D.cpp b/src/Float32ProjectionData3D.cpp index 2bd0447..680ad55 100644 --- a/src/Float32ProjectionData3D.cpp +++ b/src/Float32ProjectionData3D.cpp @@ -53,13 +53,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const CFloat32Pro CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); #ifdef _DEBUG CProjectionGeometry3D * pDataGeometry = _data.getGeometry(); - int iThisProjectionDetectorCount = pThisGeometry->getDetectorRowCount() * pThisGeometry->getDetectorColCount(); - int iDataProjectionDetectorCount = pDataGeometry->getDetectorRowCount() * pDataGeometry->getDetectorColCount(); + int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount(); ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount()); - ASTRA_ASSERT(iThisProjectionDetectorCount == iDataProjectionDetectorCount); + ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount); #endif for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) @@ -67,7 +67,7 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const CFloat32Pro CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex]; @@ -91,13 +91,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const CFloat32Pro CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); #ifdef _DEBUG CProjectionGeometry3D * pDataGeometry = _data.getGeometry(); - int iThisProjectionDetectorCount = pThisGeometry->getDetectorRowCount() * pThisGeometry->getDetectorColCount(); - int iDataProjectionDetectorCount = pDataGeometry->getDetectorRowCount() * pDataGeometry->getDetectorColCount(); + int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount(); ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount()); - ASTRA_ASSERT(iThisProjectionDetectorCount == iDataProjectionDetectorCount); + ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount); #endif for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) @@ -105,7 +105,7 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const CFloat32Pro CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex]; @@ -129,13 +129,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const CFloat32Pro CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); #ifdef _DEBUG CProjectionGeometry3D * pDataGeometry = _data.getGeometry(); - int iThisProjectionDetectorCount = pThisGeometry->getDetectorRowCount() * pThisGeometry->getDetectorColCount(); - int iDataProjectionDetectorCount = pDataGeometry->getDetectorRowCount() * pDataGeometry->getDetectorColCount(); + int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount(); ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount()); - ASTRA_ASSERT(iThisProjectionDetectorCount == iDataProjectionDetectorCount); + ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount); #endif for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) @@ -143,7 +143,7 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const CFloat32Pro CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex]; @@ -167,12 +167,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const float32& _f CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) { CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; @@ -194,12 +195,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator/=(const float32& _f CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) { CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; @@ -221,12 +223,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const float32& _f CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) { CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; @@ -248,12 +251,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const float32& _f CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) { CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; diff --git a/src/GeometryUtil3D.cpp b/src/GeometryUtil3D.cpp index 52dd5a9..c6bfd8b 100644 --- a/src/GeometryUtil3D.cpp +++ b/src/GeometryUtil3D.cpp @@ -28,8 +28,96 @@ $Id$ #include "astra/GeometryUtil3D.h" +#include <cmath> + namespace astra { + +SConeProjection* genConeProjections(unsigned int iProjAngles, + unsigned int iProjU, + unsigned int iProjV, + double fOriginSourceDistance, + double fOriginDetectorDistance, + double fDetUSize, + double fDetVSize, + const float *pfAngles) +{ + SConeProjection base; + base.fSrcX = 0.0f; + base.fSrcY = -fOriginSourceDistance; + base.fSrcZ = 0.0f; + + base.fDetSX = iProjU * fDetUSize * -0.5f; + base.fDetSY = fOriginDetectorDistance; + base.fDetSZ = iProjV * fDetVSize * -0.5f; + + base.fDetUX = fDetUSize; + base.fDetUY = 0.0f; + base.fDetUZ = 0.0f; + + base.fDetVX = 0.0f; + base.fDetVY = 0.0f; + base.fDetVZ = fDetVSize; + + SConeProjection* p = new SConeProjection[iProjAngles]; + +#define ROTATE0(name,i,alpha) do { p[i].f##name##X = base.f##name##X * cos(alpha) - base.f##name##Y * sin(alpha); p[i].f##name##Y = base.f##name##X * sin(alpha) + base.f##name##Y * cos(alpha); p[i].f##name##Z = base.f##name##Z; } while(0) + + for (unsigned int i = 0; i < iProjAngles; ++i) { + ROTATE0(Src, i, pfAngles[i]); + ROTATE0(DetS, i, pfAngles[i]); + ROTATE0(DetU, i, pfAngles[i]); + ROTATE0(DetV, i, pfAngles[i]); + } + +#undef ROTATE0 + + return p; +} + +SPar3DProjection* genPar3DProjections(unsigned int iProjAngles, + unsigned int iProjU, + unsigned int iProjV, + double fDetUSize, + double fDetVSize, + const float *pfAngles) +{ + SPar3DProjection base; + base.fRayX = 0.0f; + base.fRayY = 1.0f; + base.fRayZ = 0.0f; + + base.fDetSX = iProjU * fDetUSize * -0.5f; + base.fDetSY = 0.0f; + base.fDetSZ = iProjV * fDetVSize * -0.5f; + + base.fDetUX = fDetUSize; + base.fDetUY = 0.0f; + base.fDetUZ = 0.0f; + + base.fDetVX = 0.0f; + base.fDetVY = 0.0f; + base.fDetVZ = fDetVSize; + + SPar3DProjection* p = new SPar3DProjection[iProjAngles]; + +#define ROTATE0(name,i,alpha) do { p[i].f##name##X = base.f##name##X * cos(alpha) - base.f##name##Y * sin(alpha); p[i].f##name##Y = base.f##name##X * sin(alpha) + base.f##name##Y * cos(alpha); p[i].f##name##Z = base.f##name##Z; } while(0) + + for (unsigned int i = 0; i < iProjAngles; ++i) { + ROTATE0(Ray, i, pfAngles[i]); + ROTATE0(DetS, i, pfAngles[i]); + ROTATE0(DetU, i, pfAngles[i]); + ROTATE0(DetV, i, pfAngles[i]); + } + +#undef ROTATE0 + + return p; +} + + + + // (See declaration in header for (mathematical) description of these functions) @@ -72,4 +160,88 @@ void computeBP_UV_Coeffs(const SConeProjection& proj, double &fUX, double &fUY, } +// TODO: Handle cases of rays parallel to coordinate planes + +void backprojectPointX(const SPar3DProjection& proj, double fU, double fV, + double fX, double &fY, double &fZ) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - px) / proj.fRayX; + + fY = py + a * proj.fRayY; + fZ = pz + a * proj.fRayZ; +} + +void backprojectPointY(const SPar3DProjection& proj, double fU, double fV, + double fY, double &fX, double &fZ) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - py) / proj.fRayY; + + fX = px + a * proj.fRayX; + fZ = pz + a * proj.fRayZ; + +} + +void backprojectPointZ(const SPar3DProjection& proj, double fU, double fV, + double fZ, double &fX, double &fY) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - pz) / proj.fRayZ; + + fX = px + a * proj.fRayX; + fY = py + a * proj.fRayY; +} + + + +void backprojectPointX(const SConeProjection& proj, double fU, double fV, + double fX, double &fY, double &fZ) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - proj.fSrcX) / (px - proj.fSrcX); + + fY = proj.fSrcY + a * (py - proj.fSrcY); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); +} + +void backprojectPointY(const SConeProjection& proj, double fU, double fV, + double fY, double &fX, double &fZ) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - proj.fSrcY) / (py - proj.fSrcY); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); +} + +void backprojectPointZ(const SConeProjection& proj, double fU, double fV, + double fZ, double &fX, double &fY) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - proj.fSrcZ) / (pz - proj.fSrcZ); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fY = proj.fSrcY + a * (py - proj.fSrcY); +} + + } diff --git a/src/Globals.cpp b/src/Globals.cpp index 813f9c9..904a459 100644 --- a/src/Globals.cpp +++ b/src/Globals.cpp @@ -28,5 +28,8 @@ $Id$ #include "astra/Globals.h" +namespace astra{ + bool running_in_matlab=false; +} // nothing to see here :) diff --git a/src/Logging.cpp b/src/Logging.cpp index 8290ca0..cd7e3f0 100644 --- a/src/Logging.cpp +++ b/src/Logging.cpp @@ -70,37 +70,65 @@ void CLogger::disable() void CLogger::debug(const char *sfile, int sline, const char *fmt, ...) { _assureIsInitialized(); - va_list ap; - va_start(ap, fmt); - if(m_bEnabledScreen) clog_debug(sfile,sline,0,fmt,ap); - if(m_bEnabledFile && m_bFileProvided) clog_debug(sfile,sline,1,fmt,ap); + va_list ap, apf; + if(m_bEnabledScreen){ + va_start(ap, fmt); + clog_debug(sfile,sline,0,fmt,ap); + va_end(ap); + } + if(m_bEnabledFile && m_bFileProvided){ + va_start(apf, fmt); + clog_debug(sfile,sline,1,fmt,apf); + va_end(apf); + } } void CLogger::info(const char *sfile, int sline, const char *fmt, ...) { _assureIsInitialized(); - va_list ap; - va_start(ap, fmt); - if(m_bEnabledScreen) clog_info(sfile,sline,0,fmt,ap); - if(m_bEnabledFile && m_bFileProvided) clog_info(sfile,sline,1,fmt,ap); + va_list ap, apf; + if(m_bEnabledScreen){ + va_start(ap, fmt); + clog_info(sfile,sline,0,fmt,ap); + va_end(ap); + } + if(m_bEnabledFile && m_bFileProvided){ + va_start(apf, fmt); + clog_info(sfile,sline,1,fmt,apf); + va_end(apf); + } } void CLogger::warn(const char *sfile, int sline, const char *fmt, ...) { _assureIsInitialized(); - va_list ap; - va_start(ap, fmt); - if(m_bEnabledScreen) clog_warn(sfile,sline,0,fmt,ap); - if(m_bEnabledFile && m_bFileProvided) clog_warn(sfile,sline,1,fmt,ap); + va_list ap, apf; + if(m_bEnabledScreen){ + va_start(ap, fmt); + clog_warn(sfile,sline,0,fmt,ap); + va_end(ap); + } + if(m_bEnabledFile && m_bFileProvided){ + va_start(apf, fmt); + clog_warn(sfile,sline,1,fmt,apf); + va_end(apf); + } } void CLogger::error(const char *sfile, int sline, const char *fmt, ...) { _assureIsInitialized(); - va_list ap; - va_start(ap, fmt); - if(m_bEnabledScreen) clog_error(sfile,sline,0,fmt,ap); - if(m_bEnabledFile && m_bFileProvided) clog_error(sfile,sline,1,fmt,ap); + va_list ap, apf; + if(m_bEnabledScreen){ + va_start(ap, fmt); + clog_error(sfile,sline,0,fmt,ap); + va_end(ap); + } + if(m_bEnabledFile && m_bFileProvided){ + va_start(apf, fmt); + clog_error(sfile,sline,1,fmt,apf); + va_end(apf); + } } void CLogger::_setLevel(int id, log_level m_eLevel) diff --git a/src/ParallelProjectionGeometry3D.cpp b/src/ParallelProjectionGeometry3D.cpp index 1c87157..7b64fd9 100644 --- a/src/ParallelProjectionGeometry3D.cpp +++ b/src/ParallelProjectionGeometry3D.cpp @@ -27,8 +27,10 @@ $Id$ */ #include "astra/ParallelProjectionGeometry3D.h" -#include <boost/lexical_cast.hpp> +#include "astra/GeometryUtil3D.h" + +#include <boost/lexical_cast.hpp> #include <cstring> using namespace std; @@ -185,9 +187,9 @@ CVector3D CParallelProjectionGeometry3D::getProjectionDirection(int _iProjection return CVector3D(fDirX, fDirY, fDirZ); } -void CParallelProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, +void CParallelProjectionGeometry3D::projectPoint(double fX, double fY, double fZ, int iAngleIndex, - float32 &fU, float32 &fV) const + double &fU, double &fV) const { ASTRA_ASSERT(iAngleIndex >= 0); ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); @@ -214,6 +216,79 @@ CParallelProjectionGeometry2D * CParallelProjectionGeometry3D::createProjectionG return pOutput; } +void CParallelProjectionGeometry3D::backprojectPointX(int iAngleIndex, double fU, double fV, + double fX, double &fY, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SPar3DProjection *projs = genPar3DProjections(1, m_iDetectorColCount, m_iDetectorRowCount, + m_fDetectorSpacingX, m_fDetectorSpacingY, + &m_pfProjectionAngles[iAngleIndex]); + + SPar3DProjection &proj = projs[0]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - px) / proj.fRayX; + + fY = py + a * proj.fRayY; + fZ = pz + a * proj.fRayZ; + + delete[] projs; +} + +void CParallelProjectionGeometry3D::backprojectPointY(int iAngleIndex, double fU, double fV, + double fY, double &fX, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SPar3DProjection *projs = genPar3DProjections(1, m_iDetectorColCount, m_iDetectorRowCount, + m_fDetectorSpacingX, m_fDetectorSpacingY, + &m_pfProjectionAngles[iAngleIndex]); + + SPar3DProjection &proj = projs[0]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - py) / proj.fRayY; + + fX = px + a * proj.fRayX; + fZ = pz + a * proj.fRayZ; + + delete[] projs; +} + +void CParallelProjectionGeometry3D::backprojectPointZ(int iAngleIndex, double fU, double fV, + double fZ, double &fX, double &fY) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SPar3DProjection *projs = genPar3DProjections(1, m_iDetectorColCount, m_iDetectorRowCount, + m_fDetectorSpacingX, m_fDetectorSpacingY, + &m_pfProjectionAngles[iAngleIndex]); + + SPar3DProjection &proj = projs[0]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - pz) / proj.fRayZ; + + fX = px + a * proj.fRayX; + fY = py + a * proj.fRayY; + + delete[] projs; +} + + //---------------------------------------------------------------------------------------- } // end namespace astra diff --git a/src/ParallelVecProjectionGeometry3D.cpp b/src/ParallelVecProjectionGeometry3D.cpp index ffad6d0..d04400b 100644 --- a/src/ParallelVecProjectionGeometry3D.cpp +++ b/src/ParallelVecProjectionGeometry3D.cpp @@ -239,9 +239,9 @@ CVector3D CParallelVecProjectionGeometry3D::getProjectionDirection(int _iProject return CVector3D(p.fRayX, p.fRayY, p.fRayZ); } -void CParallelVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, - int iAngleIndex, - float32 &fU, float32 &fV) const +void CParallelVecProjectionGeometry3D::projectPoint(double fX, double fY, double fZ, + int iAngleIndex, + double &fU, double &fV) const { ASTRA_ASSERT(iAngleIndex >= 0); ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); @@ -258,6 +258,61 @@ void CParallelVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, floa } +void CParallelVecProjectionGeometry3D::backprojectPointX(int iAngleIndex, double fU, double fV, + double fX, double &fY, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SPar3DProjection &proj = m_pProjectionAngles[iAngleIndex]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - px) / proj.fRayX; + + fY = py + a * proj.fRayY; + fZ = pz + a * proj.fRayZ; +} + +void CParallelVecProjectionGeometry3D::backprojectPointY(int iAngleIndex, double fU, double fV, + double fY, double &fX, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SPar3DProjection &proj = m_pProjectionAngles[iAngleIndex]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - py) / proj.fRayY; + + fX = px + a * proj.fRayX; + fZ = pz + a * proj.fRayZ; +} + +void CParallelVecProjectionGeometry3D::backprojectPointZ(int iAngleIndex, double fU, double fV, + double fZ, double &fX, double &fY) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SPar3DProjection &proj = m_pProjectionAngles[iAngleIndex]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - pz) / proj.fRayZ; + + fX = px + a * proj.fRayX; + fY = py + a * proj.fRayY; +} + + //---------------------------------------------------------------------------------------- bool CParallelVecProjectionGeometry3D::_check() diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp new file mode 100644 index 0000000..8f7dfc5 --- /dev/null +++ b/src/PluginAlgorithm.cpp @@ -0,0 +1,401 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp + 2014-2015, CWI, Amsterdam + +Contact: astra@uantwerpen.be +Website: http://sf.net/projects/astra-toolbox + +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 <http://www.gnu.org/licenses/>. + +----------------------------------------------------------------------- +$Id$ +*/ + +#ifdef ASTRA_PYTHON + +#include "astra/PluginAlgorithm.h" +#include "astra/Logging.h" +#include <boost/algorithm/string.hpp> +#include <boost/algorithm/string/split.hpp> +#include <boost/lexical_cast.hpp> +#include <iostream> +#include <fstream> +#include <string> + +#include <Python.h> +#include "bytesobject.h" + +namespace astra { + + + +void logPythonError(){ + if(PyErr_Occurred()){ + PyObject *ptype, *pvalue, *ptraceback; + PyErr_Fetch(&ptype, &pvalue, &ptraceback); + PyErr_NormalizeException(&ptype, &pvalue, &ptraceback); + PyObject *traceback = PyImport_ImportModule("traceback"); + if(traceback!=NULL){ + PyObject *exc; + if(ptraceback==NULL){ + exc = PyObject_CallMethod(traceback,"format_exception_only","OO",ptype, pvalue); + }else{ + exc = PyObject_CallMethod(traceback,"format_exception","OOO",ptype, pvalue, ptraceback); + } + if(exc!=NULL){ + PyObject *six = PyImport_ImportModule("six"); + if(six!=NULL){ + PyObject *iter = PyObject_GetIter(exc); + if(iter!=NULL){ + PyObject *line; + std::string errStr = ""; + while(line = PyIter_Next(iter)){ + PyObject *retb = PyObject_CallMethod(six,"b","O",line); + if(retb!=NULL){ + errStr += std::string(PyBytes_AsString(retb)); + Py_DECREF(retb); + } + Py_DECREF(line); + } + ASTRA_ERROR("%s",errStr.c_str()); + Py_DECREF(iter); + } + Py_DECREF(six); + } + Py_DECREF(exc); + } + Py_DECREF(traceback); + } + if(ptype!=NULL) Py_DECREF(ptype); + if(pvalue!=NULL) Py_DECREF(pvalue); + if(ptraceback!=NULL) Py_DECREF(ptraceback); + } +} + + +CPluginAlgorithm::CPluginAlgorithm(PyObject* pyclass){ + instance = PyObject_CallObject(pyclass, NULL); + if(instance==NULL) logPythonError(); +} + +CPluginAlgorithm::~CPluginAlgorithm(){ + if(instance!=NULL){ + Py_DECREF(instance); + instance = NULL; + } +} + +bool CPluginAlgorithm::initialize(const Config& _cfg){ + if(instance==NULL) return false; + PyObject *cfgDict = XMLNode2dict(_cfg.self); + PyObject *retVal = PyObject_CallMethod(instance, "astra_init", "O",cfgDict); + Py_DECREF(cfgDict); + if(retVal==NULL){ + logPythonError(); + return false; + } + m_bIsInitialized = true; + Py_DECREF(retVal); + return m_bIsInitialized; +} + +void CPluginAlgorithm::run(int _iNrIterations){ + if(instance==NULL) return; + PyGILState_STATE state = PyGILState_Ensure(); + PyObject *retVal = PyObject_CallMethod(instance, "run", "i",_iNrIterations); + if(retVal==NULL){ + logPythonError(); + }else{ + Py_DECREF(retVal); + } + PyGILState_Release(state); +} + +void fixLapackLoading(){ + // When running in Matlab, we need to force numpy + // to use its internal lapack library instead of + // Matlab's MKL library to avoid errors. To do this, + // we set Python's dlopen flags to RTLD_NOW|RTLD_DEEPBIND + // and import 'numpy.linalg.lapack_lite' here. We reset + // Python's dlopen flags afterwards. + PyObject *sys = PyImport_ImportModule("sys"); + if(sys!=NULL){ + PyObject *curFlags = PyObject_CallMethod(sys,"getdlopenflags",NULL); + if(curFlags!=NULL){ + PyObject *retVal = PyObject_CallMethod(sys, "setdlopenflags", "i",10); + if(retVal!=NULL){ + PyObject *lapack = PyImport_ImportModule("numpy.linalg.lapack_lite"); + if(lapack!=NULL){ + Py_DECREF(lapack); + } + PyObject_CallMethod(sys, "setdlopenflags", "O",curFlags); + Py_DECREF(retVal); + } + Py_DECREF(curFlags); + } + Py_DECREF(sys); + } +} + +CPluginAlgorithmFactory::CPluginAlgorithmFactory(){ + if(!Py_IsInitialized()){ + Py_Initialize(); + PyEval_InitThreads(); + } +#ifndef _MSC_VER + if(astra::running_in_matlab) fixLapackLoading(); +#endif + pluginDict = PyDict_New(); + inspect = PyImport_ImportModule("inspect"); + six = PyImport_ImportModule("six"); +} + +CPluginAlgorithmFactory::~CPluginAlgorithmFactory(){ + if(pluginDict!=NULL){ + Py_DECREF(pluginDict); + } + if(inspect!=NULL) Py_DECREF(inspect); + if(six!=NULL) Py_DECREF(six); +} + +PyObject * getClassFromString(std::string str){ + std::vector<std::string> items; + boost::split(items, str, boost::is_any_of(".")); + PyObject *pyclass = PyImport_ImportModule(items[0].c_str()); + if(pyclass==NULL){ + logPythonError(); + return NULL; + } + PyObject *submod = pyclass; + for(unsigned int i=1;i<items.size();i++){ + submod = PyObject_GetAttrString(submod,items[i].c_str()); + Py_DECREF(pyclass); + pyclass = submod; + if(pyclass==NULL){ + logPythonError(); + return NULL; + } + } + return pyclass; +} + +bool CPluginAlgorithmFactory::registerPlugin(std::string name, std::string className){ + PyObject *str = PyBytes_FromString(className.c_str()); + PyDict_SetItemString(pluginDict, name.c_str(), str); + Py_DECREF(str); + return true; +} + +bool CPluginAlgorithmFactory::registerPlugin(std::string className){ + PyObject *pyclass = getClassFromString(className); + if(pyclass==NULL) return false; + bool ret = registerPluginClass(pyclass); + Py_DECREF(pyclass); + return ret; +} + +bool CPluginAlgorithmFactory::registerPluginClass(std::string name, PyObject * className){ + PyDict_SetItemString(pluginDict, name.c_str(), className); + return true; +} + +bool CPluginAlgorithmFactory::registerPluginClass(PyObject * className){ + PyObject *astra_name = PyObject_GetAttrString(className,"astra_name"); + if(astra_name==NULL){ + logPythonError(); + return false; + } + PyObject *retb = PyObject_CallMethod(six,"b","O",astra_name); + if(retb!=NULL){ + PyDict_SetItemString(pluginDict,PyBytes_AsString(retb),className); + Py_DECREF(retb); + }else{ + logPythonError(); + } + Py_DECREF(astra_name); + return true; +} + +CPluginAlgorithm * CPluginAlgorithmFactory::getPlugin(std::string name){ + PyObject *className = PyDict_GetItemString(pluginDict, name.c_str()); + if(className==NULL) return NULL; + CPluginAlgorithm *alg = NULL; + if(PyBytes_Check(className)){ + std::string str = std::string(PyBytes_AsString(className)); + PyObject *pyclass = getClassFromString(str); + if(pyclass!=NULL){ + alg = new CPluginAlgorithm(pyclass); + Py_DECREF(pyclass); + } + }else{ + alg = new CPluginAlgorithm(className); + } + return alg; +} + +PyObject * CPluginAlgorithmFactory::getRegistered(){ + Py_INCREF(pluginDict); + return pluginDict; +} + +std::map<std::string, std::string> CPluginAlgorithmFactory::getRegisteredMap(){ + std::map<std::string, std::string> ret; + PyObject *key, *value; + Py_ssize_t pos = 0; + while (PyDict_Next(pluginDict, &pos, &key, &value)) { + PyObject *keystr = PyObject_Str(key); + if(keystr!=NULL){ + PyObject *valstr = PyObject_Str(value); + if(valstr!=NULL){ + PyObject * keyb = PyObject_CallMethod(six,"b","O",keystr); + if(keyb!=NULL){ + PyObject * valb = PyObject_CallMethod(six,"b","O",valstr); + if(valb!=NULL){ + ret[PyBytes_AsString(keyb)] = PyBytes_AsString(valb); + Py_DECREF(valb); + } + Py_DECREF(keyb); + } + Py_DECREF(valstr); + } + Py_DECREF(keystr); + } + logPythonError(); + } + return ret; +} + +std::string CPluginAlgorithmFactory::getHelp(std::string name){ + PyObject *className = PyDict_GetItemString(pluginDict, name.c_str()); + if(className==NULL){ + ASTRA_ERROR("Plugin %s not found!",name.c_str()); + PyErr_Clear(); + return ""; + } + std::string ret = ""; + PyObject *pyclass; + if(PyBytes_Check(className)){ + std::string str = std::string(PyBytes_AsString(className)); + pyclass = getClassFromString(str); + }else{ + pyclass = className; + } + if(pyclass==NULL) return ""; + if(inspect!=NULL && six!=NULL){ + PyObject *retVal = PyObject_CallMethod(inspect,"getdoc","O",pyclass); + if(retVal!=NULL){ + if(retVal!=Py_None){ + PyObject *retb = PyObject_CallMethod(six,"b","O",retVal); + if(retb!=NULL){ + ret = std::string(PyBytes_AsString(retb)); + Py_DECREF(retb); + } + } + Py_DECREF(retVal); + }else{ + logPythonError(); + } + } + if(PyBytes_Check(className)){ + Py_DECREF(pyclass); + } + return ret; +} + +DEFINE_SINGLETON(CPluginAlgorithmFactory); + +#if PY_MAJOR_VERSION >= 3 +PyObject * pyStringFromString(std::string str){ + return PyUnicode_FromString(str.c_str()); +} +#else +PyObject * pyStringFromString(std::string str){ + return PyBytes_FromString(str.c_str()); +} +#endif + +PyObject* stringToPythonValue(std::string str){ + if(str.find(";")!=std::string::npos){ + std::vector<std::string> rows, row; + boost::split(rows, str, boost::is_any_of(";")); + PyObject *mat = PyList_New(rows.size()); + for(unsigned int i=0; i<rows.size(); i++){ + boost::split(row, rows[i], boost::is_any_of(",")); + PyObject *rowlist = PyList_New(row.size()); + for(unsigned int j=0;j<row.size();j++){ + PyList_SetItem(rowlist, j, PyFloat_FromDouble(boost::lexical_cast<double>(row[j]))); + } + PyList_SetItem(mat, i, rowlist); + } + return mat; + } + if(str.find(",")!=std::string::npos){ + std::vector<std::string> vec; + boost::split(vec, str, boost::is_any_of(",")); + PyObject *veclist = PyList_New(vec.size()); + for(unsigned int i=0;i<vec.size();i++){ + PyList_SetItem(veclist, i, PyFloat_FromDouble(boost::lexical_cast<double>(vec[i]))); + } + return veclist; + } + try{ + return PyLong_FromLong(boost::lexical_cast<long>(str)); + }catch(const boost::bad_lexical_cast &){ + try{ + return PyFloat_FromDouble(boost::lexical_cast<double>(str)); + }catch(const boost::bad_lexical_cast &){ + return pyStringFromString(str); + } + } +} + +PyObject* XMLNode2dict(XMLNode node){ + PyObject *dct = PyDict_New(); + PyObject *opts = PyDict_New(); + if(node.hasAttribute("type")){ + PyObject *obj = pyStringFromString(node.getAttribute("type").c_str()); + PyDict_SetItemString(dct, "type", obj); + Py_DECREF(obj); + } + std::list<XMLNode> nodes = node.getNodes(); + std::list<XMLNode>::iterator it = nodes.begin(); + while(it!=nodes.end()){ + XMLNode subnode = *it; + if(subnode.getName()=="Option"){ + PyObject *obj; + if(subnode.hasAttribute("value")){ + obj = stringToPythonValue(subnode.getAttribute("value")); + }else{ + obj = stringToPythonValue(subnode.getContent()); + } + PyDict_SetItemString(opts, subnode.getAttribute("key").c_str(), obj); + Py_DECREF(obj); + }else{ + PyObject *obj = stringToPythonValue(subnode.getContent()); + PyDict_SetItemString(dct, subnode.getName().c_str(), obj); + Py_DECREF(obj); + } + ++it; + } + PyDict_SetItemString(dct, "options", opts); + Py_DECREF(opts); + return dct; +} + +} +#endif diff --git a/src/ReconstructionAlgorithm2D.cpp b/src/ReconstructionAlgorithm2D.cpp index 767efe6..4575ff7 100644 --- a/src/ReconstructionAlgorithm2D.cpp +++ b/src/ReconstructionAlgorithm2D.cpp @@ -85,9 +85,16 @@ bool CReconstructionAlgorithm2D::initialize(const Config& _cfg) // projector XMLNode node = _cfg.self.getSingleNode("ProjectorId"); - ASTRA_CONFIG_CHECK(node, "Reconstruction2D", "No ProjectorId tag specified."); - int id = boost::lexical_cast<int>(node.getContent()); - m_pProjector = CProjector2DManager::getSingleton().get(id); + if (requiresProjector()) { + ASTRA_CONFIG_CHECK(node, "Reconstruction2D", "No ProjectorId tag specified."); + } + int id; + if (node) { + id = boost::lexical_cast<int>(node.getContent()); + m_pProjector = CProjector2DManager::getSingleton().get(id); + } else { + m_pProjector = 0; + } CC.markNodeParsed("ProjectorId"); // sinogram data @@ -205,18 +212,22 @@ void CReconstructionAlgorithm2D::setSinogramMask(CFloat32ProjectionData2D* _pMas bool CReconstructionAlgorithm2D::_check() { // check pointers - ASTRA_CONFIG_CHECK(m_pProjector, "Reconstruction2D", "Invalid Projector Object."); + if (requiresProjector()) + ASTRA_CONFIG_CHECK(m_pProjector, "Reconstruction2D", "Invalid Projector Object."); ASTRA_CONFIG_CHECK(m_pSinogram, "Reconstruction2D", "Invalid Projection Data Object."); ASTRA_CONFIG_CHECK(m_pReconstruction, "Reconstruction2D", "Invalid Reconstruction Data Object."); // check initializations - ASTRA_CONFIG_CHECK(m_pProjector->isInitialized(), "Reconstruction2D", "Projector Object Not Initialized."); + if (requiresProjector()) + ASTRA_CONFIG_CHECK(m_pProjector->isInitialized(), "Reconstruction2D", "Projector Object Not Initialized."); ASTRA_CONFIG_CHECK(m_pSinogram->isInitialized(), "Reconstruction2D", "Projection Data Object Not Initialized."); ASTRA_CONFIG_CHECK(m_pReconstruction->isInitialized(), "Reconstruction2D", "Reconstruction Data Object Not Initialized."); // check compatibility between projector and data classes - ASTRA_CONFIG_CHECK(m_pSinogram->getGeometry()->isEqual(m_pProjector->getProjectionGeometry()), "Reconstruction2D", "Projection Data not compatible with the specified Projector."); - ASTRA_CONFIG_CHECK(m_pReconstruction->getGeometry()->isEqual(m_pProjector->getVolumeGeometry()), "Reconstruction2D", "Reconstruction Data not compatible with the specified Projector."); + if (requiresProjector()) { + ASTRA_CONFIG_CHECK(m_pSinogram->getGeometry()->isEqual(m_pProjector->getProjectionGeometry()), "Reconstruction2D", "Projection Data not compatible with the specified Projector."); + ASTRA_CONFIG_CHECK(m_pReconstruction->getGeometry()->isEqual(m_pProjector->getVolumeGeometry()), "Reconstruction2D", "Reconstruction Data not compatible with the specified Projector."); + } // success return true; diff --git a/src/VolumeGeometry3D.cpp b/src/VolumeGeometry3D.cpp index a1cf424..3de146f 100644 --- a/src/VolumeGeometry3D.cpp +++ b/src/VolumeGeometry3D.cpp @@ -45,6 +45,7 @@ bool CVolumeGeometry3D::_check() ASTRA_CONFIG_CHECK(m_fWindowMinZ < m_fWindowMaxZ, "VolumeGeometry3D", "WindowMinZ should be lower than WindowMaxZ."); ASTRA_CONFIG_CHECK(m_iGridTotCount == (m_iGridColCount * m_iGridRowCount * m_iGridSliceCount), "VolumeGeometry3D", "Internal configuration error."); +#if 0 ASTRA_CONFIG_CHECK(m_fWindowLengthX == (m_fWindowMaxX - m_fWindowMinX), "VolumeGeometry3D", "Internal configuration error."); ASTRA_CONFIG_CHECK(m_fWindowLengthY == (m_fWindowMaxY - m_fWindowMinY), "VolumeGeometry3D", "Internal configuration error."); ASTRA_CONFIG_CHECK(m_fWindowLengthZ == (m_fWindowMaxZ - m_fWindowMinZ), "VolumeGeometry3D", "Internal configuration error."); @@ -57,6 +58,7 @@ bool CVolumeGeometry3D::_check() ASTRA_CONFIG_CHECK(m_fDivPixelLengthX == (1.0f / m_fPixelLengthX), "VolumeGeometry3D", "Internal configuration error."); ASTRA_CONFIG_CHECK(m_fDivPixelLengthY == (1.0f / m_fPixelLengthY), "VolumeGeometry3D", "Internal configuration error."); ASTRA_CONFIG_CHECK(m_fDivPixelLengthZ == (1.0f / m_fPixelLengthZ), "VolumeGeometry3D", "Internal configuration error."); +#endif return true; } |