Add CUDA parvec support

7 files changed, 242 insertions, 326 deletions

diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp
index aa97eec..9d23253 100644
--- a/src/CudaFilteredBackProjectionAlgorithm.cpp
+++ b/src/CudaFilteredBackProjectionAlgorithm.cpp
@@ -33,6 +33,7 @@ $Id$
 #include "astra/AstraObjectManager.h"
 #include "astra/CudaProjector2D.h"
 #include "../cuda/2d/astra.h"
+#include "../cuda/2d/fbp.h"
 
 #include "astra/Logging.h"
 
@@ -44,8 +45,7 @@ string CCudaFilteredBackProjectionAlgorithm::type = "FBP_CUDA";
 CCudaFilteredBackProjectionAlgorithm::CCudaFilteredBackProjectionAlgorithm()
 {
 	m_bIsInitialized = false;
-	CReconstructionAlgorithm2D::_clear();
-	m_pFBP = 0;
+	CCudaReconstructionAlgorithm2D::_clear();
 	m_pfFilter = NULL;
 	m_fFilterParameter = -1.0f;
 	m_fFilterD = 1.0f;
@@ -53,35 +53,8 @@ CCudaFilteredBackProjectionAlgorithm::CCudaFilteredBackProjectionAlgorithm()
 
 CCudaFilteredBackProjectionAlgorithm::~CCudaFilteredBackProjectionAlgorithm()
 {
-	if(m_pfFilter != NULL)
-	{
-		delete [] m_pfFilter;
-		m_pfFilter = NULL;
-	}
-
-	if(m_pFBP != NULL)
-	{
-		delete m_pFBP;
-		m_pFBP = NULL;
-	}
-}
-
-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();
-	}
-
+	delete[] m_pfFilter;
+	m_pfFilter = NULL;
 }
 
 bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
@@ -92,54 +65,28 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
 	// if already initialized, clear first
 	if (m_bIsInitialized)
 	{
+#warning FIXME Necessary?
 		clear();
 	}
 
-	// Projector
-	XMLNode node = _cfg.self.getSingleNode("ProjectorId");
-	CCudaProjector2D* pCudaProjector = 0;
-	if (node) {
-		int id = node.getContentInt();
-		CProjector2D *projector = CProjector2DManager::getSingleton().get(id);
-		pCudaProjector = dynamic_cast<CCudaProjector2D*>(projector);
-		if (!pCudaProjector) {
-			ASTRA_WARN("non-CUDA Projector2D passed");
-		}
-	}
-	CC.markNodeParsed("ProjectorId");
-
-
-	// sinogram data
-	node = _cfg.self.getSingleNode("ProjectionDataId");
-	ASTRA_CONFIG_CHECK(node, "CudaFBP", "No ProjectionDataId tag specified.");
-	int id = node.getContentInt();
-	m_pSinogram = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id));
-	CC.markNodeParsed("ProjectionDataId");
+	m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_cfg);
+	if (!m_bIsInitialized)
+		return false;
 
-	// reconstruction data
-	node = _cfg.self.getSingleNode("ReconstructionDataId");
-	ASTRA_CONFIG_CHECK(node, "CudaFBP", "No ReconstructionDataId tag specified.");
-	id = node.getContentInt();
-	m_pReconstruction = dynamic_cast<CFloat32VolumeData2D*>(CData2DManager::getSingleton().get(id));
-	CC.markNodeParsed("ReconstructionDataId");
 
 	// filter type
-	node = _cfg.self.getSingleNode("FilterType");
+	XMLNode node = _cfg.self.getSingleNode("FilterType");
 	if (node)
-	{
 		m_eFilter = _convertStringToFilter(node.getContent().c_str());
-	}
 	else
-	{
 		m_eFilter = FILTER_RAMLAK;
-	}
 	CC.markNodeParsed("FilterType");
 
 	// filter
 	node = _cfg.self.getSingleNode("FilterSinogramId");
 	if (node)
 	{
-		id = node.getContentInt();
+		int id = node.getContentInt();
 		const CFloat32ProjectionData2D * pFilterData = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id));
 		m_iFilterWidth = pFilterData->getGeometry()->getDetectorCount();
 		int iFilterProjectionCount = pFilterData->getGeometry()->getProjectionAngleCount();
@@ -188,20 +135,9 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
 
 	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;
+	m_pAlgo = new astraCUDA::FBP();
+	m_bAlgoInit = false;
 
 	return check();
 }
@@ -226,9 +162,8 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D *
 	// success
 	m_bIsInitialized = true;
 
-	m_pFBP = new AstraFBP;
-
-	m_bAstraFBPInit = false;
+	m_pAlgo = new astraCUDA::FBP();
+	m_bAlgoInit = false;
 
 	if(_pfFilter != NULL)
 	{
@@ -256,107 +191,26 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D *
 	return check();
 }
 
-void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */)
-{
-	// check initialized
-	ASTRA_ASSERT(m_bIsInitialized);
-
-	if (!m_bAstraFBPInit) {
-
-		const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry();
-		const CParallelProjectionGeometry2D* parprojgeom = dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry());
-		const CFanFlatProjectionGeometry2D* fanprojgeom = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pSinogram->getGeometry());
-
-		bool ok = true;
-
-		// TODO: non-square pixels?
-		ok &= m_pFBP->setReconstructionGeometry(volgeom.getGridColCount(),
-		                                         volgeom.getGridRowCount(),
-		                                         volgeom.getPixelLengthX());
-		int iDetectorCount;
-		if (parprojgeom) {
-
-			float *offsets, *angles, detSize, outputScale;
-
-			ok = convertAstraGeometry(&volgeom, parprojgeom, offsets, angles, detSize, outputScale);
-
-
-			ok &= m_pFBP->setProjectionGeometry(parprojgeom->getProjectionAngleCount(),
-			                                     parprojgeom->getDetectorCount(),
-			                                     angles,
-			                                     parprojgeom->getDetectorWidth());
-			iDetectorCount = parprojgeom->getDetectorCount();
 
-			// TODO: Are detSize and outputScale handled correctly?
-
-			if (offsets)
-				ok &= m_pFBP->setTOffsets(offsets);
-			ASTRA_ASSERT(ok);
-
-			delete[] offsets;
-			delete[] angles;
-
-		} else if (fanprojgeom) {
-
-			astraCUDA::SFanProjection* projs;
-			float outputScale;
-
-			// FIXME: Implement this, and clean up the interface to AstraFBP.
-			if (abs(volgeom.getWindowMinX() + volgeom.getWindowMaxX()) > 0.00001 * volgeom.getPixelLengthX()) {
-				// Off-center volume geometry isn't supported yet
-				ASTRA_ASSERT(false);
-			}
-			if (abs(volgeom.getWindowMinY() + volgeom.getWindowMaxY()) > 0.00001 * volgeom.getPixelLengthY()) {
-				// Off-center volume geometry isn't supported yet
-				ASTRA_ASSERT(false);
-			}
-
-			ok = convertAstraGeometry(&volgeom, fanprojgeom, projs, outputScale);
-
-			// CHECKME: outputScale?
-
-			ok &= m_pFBP->setFanGeometry(fanprojgeom->getProjectionAngleCount(),
-			                             fanprojgeom->getDetectorCount(),
-			                             projs,
-			                             fanprojgeom->getProjectionAngles(),
-			                             fanprojgeom->getOriginSourceDistance(),
-			                             fanprojgeom->getOriginDetectorDistance(),
-
-		                                     fanprojgeom->getDetectorWidth(),
-			                             m_bShortScan);
-
-			iDetectorCount = fanprojgeom->getDetectorCount();
-
-			delete[] projs;
-		} else {
-			assert(false);
-		}
-
-		ok &= m_pFBP->setPixelSuperSampling(m_iPixelSuperSampling);
-
-		ASTRA_ASSERT(ok);
+void CCudaFilteredBackProjectionAlgorithm::initCUDAAlgorithm()
+{
+	CCudaReconstructionAlgorithm2D::initCUDAAlgorithm();
 
-		ok &= m_pFBP->init(m_iGPUIndex);
-		ASTRA_ASSERT(ok);
+	astraCUDA::FBP* pFBP = dynamic_cast<astraCUDA::FBP*>(m_pAlgo);
 
-		ok &= m_pFBP->setSinogram(m_pSinogram->getDataConst(), iDetectorCount);
+	bool ok = pFBP->setFilter(m_eFilter, m_pfFilter, m_iFilterWidth, m_fFilterD, m_fFilterParameter);
+	if (!ok) {
+		ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set filter");
 		ASTRA_ASSERT(ok);
-
-		ok &= m_pFBP->setFilter(m_eFilter, m_pfFilter, m_iFilterWidth, m_fFilterD, m_fFilterParameter);
-		ASTRA_ASSERT(ok);
-
-		m_bAstraFBPInit = true;
 	}
 
-	bool ok = m_pFBP->run();
-	ASTRA_ASSERT(ok);
-
-	const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry();
-	ok &= m_pFBP->getReconstruction(m_pReconstruction->getData(), volgeom.getGridColCount());
-
-	ASTRA_ASSERT(ok);
+	ok &= pFBP->setShortScan(m_bShortScan);
+	if (!ok) {
+		ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set short-scan mode");
+	}
 }
 
+
 bool CCudaFilteredBackProjectionAlgorithm::check()
 {
 	// check pointers
@@ -395,16 +249,6 @@ static int calcNextPowerOfTwo(int _iValue)
 	return iOutput;
 }
 
-int CCudaFilteredBackProjectionAlgorithm::calcIdealRealFilterWidth(int _iDetectorCount)
-{
-	return calcNextPowerOfTwo(_iDetectorCount);
-}
-
-int CCudaFilteredBackProjectionAlgorithm::calcIdealFourierFilterWidth(int _iDetectorCount)
-{
-	return (calcNextPowerOfTwo(_iDetectorCount) / 2 + 1);
-}
-
 static bool stringCompareLowerCase(const char * _stringA, const char * _stringB)
 {
 	int iCmpReturn = 0;
@@ -518,12 +362,3 @@ E_FBPFILTER CCudaFilteredBackProjectionAlgorithm::_convertStringToFilter(const c
 	return output;
 }
 
-void CCudaFilteredBackProjectionAlgorithm::testGenFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount)
-{
-	genFilter(_eFilter, _fD, _iProjectionCount, _pFilter, _iFFTRealDetectorCount, _iFFTFourierDetectorCount);
-}
-
-int CCudaFilteredBackProjectionAlgorithm::getGPUCount()
-{
-	return 0;
-}
diff --git a/src/CudaForwardProjectionAlgorithm.cpp b/src/CudaForwardProjectionAlgorithm.cpp
index 80f2e02..bdd7dd0 100644
--- a/src/CudaForwardProjectionAlgorithm.cpp
+++ b/src/CudaForwardProjectionAlgorithm.cpp
@@ -221,56 +221,48 @@ void CCudaForwardProjectionAlgorithm::run(int)
 	// check initialized
 	assert(m_bIsInitialized);
 
-	CVolumeGeometry2D* pVolGeom = m_pVolume->getGeometry();
-	const CParallelProjectionGeometry2D* parProjGeom = dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry());
-	const CFanFlatProjectionGeometry2D* fanProjGeom = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pSinogram->getGeometry());
-	const CFanFlatVecProjectionGeometry2D* fanVecProjGeom = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pSinogram->getGeometry());
+	bool ok;
 
-	bool ok = false;
-	if (parProjGeom) {
+	const CVolumeGeometry2D* pVolGeom = m_pVolume->getGeometry();
+	const CProjectionGeometry2D* pProjGeom = m_pSinogram->getGeometry();
+	astraCUDA::SDimensions dims;
 
-		float *offsets, *angles, detSize, outputScale;
-		ok = convertAstraGeometry(pVolGeom, parProjGeom, offsets, angles, detSize, outputScale);
+	ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims);
 
-		ASTRA_ASSERT(ok); // FIXME
+	if (!ok)
+		return;
 
-		// FIXME: Output scaling
+	astraCUDA::SParProjection* pParProjs = 0;
+	astraCUDA::SFanProjection* pFanProjs = 0;
+	float fOutputScale = 1.0f;
 
-		ok = astraCudaFP(m_pVolume->getDataConst(), m_pSinogram->getData(),
-		                 pVolGeom->getGridColCount(), pVolGeom->getGridRowCount(),
-		                 parProjGeom->getProjectionAngleCount(),
-		                 parProjGeom->getDetectorCount(),
-		                 angles, offsets, detSize,
-		                 m_iDetectorSuperSampling, 1.0f * outputScale, m_iGPUIndex);
-
-		delete[] offsets;
-		delete[] angles;
+	ok = convertAstraGeometry(pVolGeom, pProjGeom, pParProjs, pFanProjs, fOutputScale);
+	if (!ok)
+		return;
 
-	} else if (fanProjGeom || fanVecProjGeom) {
+	if (pParProjs) {
+		assert(!pFanProjs);
 
-		astraCUDA::SFanProjection* projs;
-		float outputScale;
+		ok = astraCudaFP(m_pVolume->getDataConst(), m_pSinogram->getData(),
+		                 pVolGeom->getGridColCount(), pVolGeom->getGridRowCount(),
+		                 pProjGeom->getProjectionAngleCount(),
+		                 pProjGeom->getDetectorCount(),
+		                 pParProjs,
+		                 m_iDetectorSuperSampling, 1.0f * fOutputScale, m_iGPUIndex);
 
-		if (fanProjGeom) {
-			ok = convertAstraGeometry(pVolGeom, fanProjGeom, projs, outputScale);
-		} else {
-			ok = convertAstraGeometry(pVolGeom, fanVecProjGeom, projs, outputScale);
-		}
+		delete[] pParProjs;
 
-		ASTRA_ASSERT(ok);
+	} else {
+		assert(pFanProjs);
 
 		ok = astraCudaFanFP(m_pVolume->getDataConst(), m_pSinogram->getData(),
 		                    pVolGeom->getGridColCount(), pVolGeom->getGridRowCount(),
-		                    m_pSinogram->getGeometry()->getProjectionAngleCount(),
-		                    m_pSinogram->getGeometry()->getDetectorCount(),
-		                    projs,
-		                    m_iDetectorSuperSampling, outputScale, m_iGPUIndex);
-
-		delete[] projs;
-
-	} else {
+		                    pProjGeom->getProjectionAngleCount(),
+		                    pProjGeom->getDetectorCount(),
+		                    pFanProjs,
+		                    m_iDetectorSuperSampling, fOutputScale, m_iGPUIndex);
 
-		ASTRA_ASSERT(false);
+		delete[] pFanProjs;
 
 	}
 
diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp
index 2798434..326ef1e 100644
--- a/src/CudaReconstructionAlgorithm2D.cpp
+++ b/src/CudaReconstructionAlgorithm2D.cpp
@@ -250,69 +250,14 @@ bool CCudaReconstructionAlgorithm2D::setupGeometry()
 	ok = m_pAlgo->setGPUIndex(m_iGPUIndex);
 	if (!ok) return false;
 
-	astraCUDA::SDimensions dims;
-
 	const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry();
+	const CProjectionGeometry2D& projgeom = *m_pSinogram->getGeometry();
 
-	// TODO: non-square pixels?
-	dims.iVolWidth = volgeom.getGridColCount();
-	dims.iVolHeight = volgeom.getGridRowCount();
-	float fPixelSize = volgeom.getPixelLengthX();
-
-	dims.iRaysPerDet = m_iDetectorSuperSampling;
-	dims.iRaysPerPixelDim = m_iPixelSuperSampling;
-
-
-	const CParallelProjectionGeometry2D* parProjGeom = dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry());
-	const CFanFlatProjectionGeometry2D* fanProjGeom = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pSinogram->getGeometry());
-	const CFanFlatVecProjectionGeometry2D* fanVecProjGeom = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pSinogram->getGeometry());
-
-	if (parProjGeom) {
-
-		float *offsets, *angles, detSize, outputScale;
-
-		ok = convertAstraGeometry(&volgeom, parProjGeom, offsets, angles, detSize, outputScale);
-
-		dims.iProjAngles = parProjGeom->getProjectionAngleCount();
-		dims.iProjDets = parProjGeom->getDetectorCount();
-		dims.fDetScale = parProjGeom->getDetectorWidth() / fPixelSize;
-
-		ok = m_pAlgo->setGeometry(dims, parProjGeom->getProjectionAngles());
-		ok &= m_pAlgo->setTOffsets(offsets);
-
-		// CHECKME: outputScale? detSize?
-
-		delete[] offsets;
-		delete[] angles;
-
-	} else if (fanProjGeom || fanVecProjGeom) {
-
-		astraCUDA::SFanProjection* projs;
-		float outputScale;
-
-		if (fanProjGeom) {
-			ok = convertAstraGeometry(&volgeom, fanProjGeom, projs, outputScale);
-		} else {
-			ok = convertAstraGeometry(&volgeom, fanVecProjGeom, projs, outputScale);
-		}
-
-		dims.iProjAngles = m_pSinogram->getGeometry()->getProjectionAngleCount();
-		dims.iProjDets = m_pSinogram->getGeometry()->getDetectorCount();
-		dims.fDetScale = m_pSinogram->getGeometry()->getDetectorWidth() / fPixelSize;
-
-		ok = m_pAlgo->setFanGeometry(dims, projs);
-
-		// CHECKME: outputScale?
-
-		delete[] projs;
-
-	} else {
-
-		ASTRA_ASSERT(false);
-
-	}
+	ok = m_pAlgo->setGeometry(&volgeom, &projgeom);
 	if (!ok) return false;
 
+	ok = m_pAlgo->setSuperSampling(m_iDetectorSuperSampling, m_iPixelSuperSampling);
+	if (!ok) return false;
 
 	if (m_bUseReconstructionMask)
 		ok &= m_pAlgo->enableVolumeMask();
diff --git a/src/FanFlatProjectionGeometry2D.cpp b/src/FanFlatProjectionGeometry2D.cpp
index 3aab582..4eec9c4 100644
--- a/src/FanFlatProjectionGeometry2D.cpp
+++ b/src/FanFlatProjectionGeometry2D.cpp
@@ -28,6 +28,8 @@ $Id$
 
 #include "astra/FanFlatProjectionGeometry2D.h"
 
+#include "astra/GeometryUtil2D.h"
+
 #include <cstring>
 #include <sstream>
 
@@ -218,16 +220,12 @@ Config* CFanFlatProjectionGeometry2D::getConfiguration() const
 //----------------------------------------------------------------------------------------
 CFanFlatVecProjectionGeometry2D* CFanFlatProjectionGeometry2D::toVectorGeometry()
 {
-	SFanProjection* vectors = new SFanProjection[m_iProjectionAngleCount];
-	for (int i = 0; i < m_iProjectionAngleCount; ++i)
-	{
-		vectors[i].fSrcX =  sinf(m_pfProjectionAngles[i]) * m_fOriginSourceDistance;
-		vectors[i].fSrcY = -cosf(m_pfProjectionAngles[i]) * m_fOriginSourceDistance;
-		vectors[i].fDetUX = cosf(m_pfProjectionAngles[i]) * m_fDetectorWidth;
-		vectors[i].fDetUY = sinf(m_pfProjectionAngles[i]) * m_fDetectorWidth;
-		vectors[i].fDetSX = -sinf(m_pfProjectionAngles[i]) * m_fOriginDetectorDistance - 0.5f * m_iDetectorCount * vectors[i].fDetUX;
-		vectors[i].fDetSY =  cosf(m_pfProjectionAngles[i]) * m_fOriginDetectorDistance - 0.5f * m_iDetectorCount * vectors[i].fDetUY;
-	}
+	SFanProjection* vectors = genFanProjections(m_iProjectionAngleCount,
+	                                            m_iDetectorCount,
+	                                            m_fOriginSourceDistance,
+	                                            m_fOriginDetectorDistance,
+	                                            m_fDetectorWidth,
+	                                            m_pfProjectionAngles);
 
 	CFanFlatVecProjectionGeometry2D* vecGeom = new CFanFlatVecProjectionGeometry2D();
 	vecGeom->initialize(m_iProjectionAngleCount, m_iDetectorCount, vectors);
diff --git a/src/GeometryUtil2D.cpp b/src/GeometryUtil2D.cpp
new file mode 100644
index 0000000..1bca2bc
--- /dev/null
+++ b/src/GeometryUtil2D.cpp
@@ -0,0 +1,161 @@
+/*
+-----------------------------------------------------------------------
+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$
+*/
+
+#include "astra/GeometryUtil2D.h"
+
+#include <cmath>
+
+namespace astra {
+
+SParProjection* genParProjections(unsigned int iProjAngles,
+                                  unsigned int iProjDets,
+                                  double fDetSize,
+                                  const float *pfAngles,
+                                  const float *pfExtraOffsets)
+{
+	SParProjection base;
+	base.fRayX = 0.0f;
+	base.fRayY = 1.0f;
+
+	base.fDetSX = iProjDets * fDetSize * -0.5f;
+	base.fDetSY = 0.0f;
+
+	base.fDetUX = fDetSize;
+	base.fDetUY = 0.0f;
+
+	SParProjection* p = new SParProjection[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); } while(0)
+
+	for (unsigned int i = 0; i < iProjAngles; ++i) {
+		if (pfExtraOffsets) {
+			// TODO
+		}
+
+		ROTATE0(Ray, i, pfAngles[i]);
+		ROTATE0(DetS, i, pfAngles[i]);
+		ROTATE0(DetU, i, pfAngles[i]);
+
+		if (pfExtraOffsets) {
+			float d = pfExtraOffsets[i];
+			p[i].fDetSX -= d * p[i].fDetUX;
+			p[i].fDetSY -= d * p[i].fDetUY;
+		}
+	}
+
+#undef ROTATE0
+
+	return p;
+}
+
+
+SFanProjection* genFanProjections(unsigned int iProjAngles,
+                                  unsigned int iProjDets,
+                                  double fOriginSource, double fOriginDetector,
+                                  double fDetSize,
+                                  const float *pfAngles)
+//                                  const float *pfExtraOffsets)
+{
+	SFanProjection *pProjs = new SFanProjection[iProjAngles];
+
+	float fSrcX0 = 0.0f;
+	float fSrcY0 = -fOriginSource;
+	float fDetUX0 = fDetSize;
+	float fDetUY0 = 0.0f;
+	float fDetSX0 = iProjDets * fDetUX0 / -2.0f;
+	float fDetSY0 = fOriginDetector;
+
+#define ROTATE0(name,i,alpha) do { pProjs[i].f##name##X = f##name##X0 * cos(alpha) - f##name##Y0 * sin(alpha); pProjs[i].f##name##Y = f##name##X0 * sin(alpha) + f##name##Y0 * cos(alpha); } 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]);
+	}
+
+#undef ROTATE0
+
+	return pProjs;
+}
+
+// Convert a SParProjection back into its set of "standard" circular parallel
+// beam parameters. This is always possible.
+bool getParParameters(const SParProjection &proj /* , ... */)
+{
+	// angle
+	// det size
+	// offset
+
+	// (see convertAndUploadAngles in par_fp.cu)
+
+	return true;
+}
+
+// Convert a SFanProjection back into its set of "standard" circular fan beam
+// parameters. This will return false if it can not be represented in this way.
+bool getFanParameters(const SFanProjection &proj, unsigned int iProjDets, float &fAngle, float &fOriginSource, float &fOriginDetector, float &fDetSize, float &fOffset)
+{
+	// angle
+	// det size
+	// offset
+	// origin-source
+	// origin-detector
+
+	// Need to check if line source-origin is orthogonal to vector ux,uy
+	// (including the case source==origin)
+
+	// (equivalent: source and origin project to same point on detector)
+
+	double dp = proj.fSrcX * proj.fDetUX + proj.fSrcY * proj.fDetUY;
+
+	double rel = (proj.fSrcX*proj.fSrcX + proj.fSrcY*proj.fSrcY) * (proj.fDetUX*proj.fDetUX + proj.fDetUY*proj.fDetUY);
+	rel = sqrt(rel);
+
+	if (std::abs(dp) > rel * 0.0001)
+		return false;
+
+	fOriginSource = sqrt(proj.fSrcX*proj.fSrcX + proj.fSrcY*proj.fSrcY);
+
+	fDetSize = sqrt(proj.fDetUX*proj.fDetUX + proj.fDetUY*proj.fDetUY);
+
+	// project origin on detector line ( == project source on detector line)
+
+	double t = (- proj.fDetSX) * proj.fDetUX + (- proj.fDetSY) * proj.fDetUY;
+
+	fOffset = (float)t - 0.5*iProjDets;
+
+	// TODO: CHECKME
+	fOriginDetector = sqrt((proj.fDetSX + t * proj.fDetUX)*(proj.fDetSX + t * proj.fDetUX) + (proj.fDetSY + t * proj.fDetUY)*(proj.fDetSY + t * proj.fDetUY));
+
+	//float fAngle = atan2(proj.fDetSX + t * proj.fDetUX - proj.fSrcX, proj.fDetSY + t * proj.fDetUY); // TODO: Fix order + sign
+	fAngle = atan2(proj.fDetUY, proj.fDetUX); // TODO: Check order + sign
+
+	return true;
+}
+
+
+}
diff --git a/src/ParallelProjectionGeometry2D.cpp b/src/ParallelProjectionGeometry2D.cpp
index 105e24b..e1f93aa 100644
--- a/src/ParallelProjectionGeometry2D.cpp
+++ b/src/ParallelProjectionGeometry2D.cpp
@@ -28,6 +28,8 @@ $Id$
 
 #include "astra/ParallelProjectionGeometry2D.h"
 
+#include "astra/GeometryUtil2D.h"
+
 #include <cstring>
 
 using namespace std;
@@ -48,15 +50,13 @@ CParallelProjectionGeometry2D::CParallelProjectionGeometry2D() :
 CParallelProjectionGeometry2D::CParallelProjectionGeometry2D(int _iProjectionAngleCount, 
 															 int _iDetectorCount, 
 															 float32 _fDetectorWidth, 
-															 const float32* _pfProjectionAngles,
-															 const float32* _pfExtraDetectorOffsets)
+															 const float32* _pfProjectionAngles)
 {
 	_clear();
 	initialize(_iProjectionAngleCount,
 				_iDetectorCount, 
 				_fDetectorWidth, 
-				_pfProjectionAngles,
-				_pfExtraDetectorOffsets);
+				_pfProjectionAngles);
 }
 
 //----------------------------------------------------------------------------------------
@@ -115,14 +115,12 @@ bool CParallelProjectionGeometry2D::initialize(const Config& _cfg)
 bool CParallelProjectionGeometry2D::initialize(int _iProjectionAngleCount, 
 											   int _iDetectorCount, 
 											   float32 _fDetectorWidth, 
-											   const float32* _pfProjectionAngles,
-											   const float32* _pfExtraDetectorOffsets)
+											   const float32* _pfProjectionAngles)
 {
 	_initialize(_iProjectionAngleCount, 
 			    _iDetectorCount, 
 			    _fDetectorWidth, 
-			    _pfProjectionAngles,
-				_pfExtraDetectorOffsets);
+			    _pfProjectionAngles);
 
 	// success
 	m_bInitialized = _check();
@@ -178,11 +176,6 @@ Config* CParallelProjectionGeometry2D::getConfiguration() const
 	cfg->self.addChildNode("DetectorCount", getDetectorCount());
 	cfg->self.addChildNode("DetectorWidth", getDetectorWidth());
 	cfg->self.addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount);
-	if(m_pfExtraDetectorOffset!=NULL){
-		XMLNode opt = cfg->self.addChildNode("Option");
-		opt.addAttribute("key","ExtraDetectorOffset");
-		opt.setContent(m_pfExtraDetectorOffset, m_iProjectionAngleCount);
-	}
 	return cfg;
 }
 
@@ -203,17 +196,11 @@ CVector3D CParallelProjectionGeometry2D::getProjectionDirection(int _iProjection
 //----------------------------------------------------------------------------------------
 CParallelVecProjectionGeometry2D* CParallelProjectionGeometry2D::toVectorGeometry()
 {
-	SParProjection* vectors = new SParProjection[m_iProjectionAngleCount];
-	for (int i = 0; i < m_iProjectionAngleCount; ++i)
-	{
-		vectors[i].fRayX = sinf(m_pfProjectionAngles[i]);
-		vectors[i].fRayY = -cosf(m_pfProjectionAngles[i]);
-		vectors[i].fDetUX = cosf(m_pfProjectionAngles[i]) * m_fDetectorWidth;
-		vectors[i].fDetUY = sinf(m_pfProjectionAngles[i]) * m_fDetectorWidth;		
-		vectors[i].fDetSX = -0.5f * m_iDetectorCount * vectors[i].fDetUX;
-		vectors[i].fDetSY = -0.5f * m_iDetectorCount * vectors[i].fDetUY;
-	}
-
+	SParProjection* vectors = genParProjections(m_iProjectionAngleCount,
+	                                            m_iDetectorCount,
+	                                            m_fDetectorWidth,
+	                                            m_pfProjectionAngles, 0);
+	// TODO: ExtraOffsets?
 	CParallelVecProjectionGeometry2D* vecGeom = new CParallelVecProjectionGeometry2D();
 	vecGeom->initialize(m_iProjectionAngleCount, m_iDetectorCount, vectors);
 	delete[] vectors;
diff --git a/src/ProjectionGeometry2D.cpp b/src/ProjectionGeometry2D.cpp
index 03fbd22..a306b48 100644
--- a/src/ProjectionGeometry2D.cpp
+++ b/src/ProjectionGeometry2D.cpp
@@ -45,11 +45,10 @@ CProjectionGeometry2D::CProjectionGeometry2D() : configCheckData(0)
 CProjectionGeometry2D::CProjectionGeometry2D(int _iAngleCount, 
 											 int _iDetectorCount, 
 											 float32 _fDetectorWidth, 
-											 const float32* _pfProjectionAngles,
-											 const float32* _pfExtraDetectorOffsets) : configCheckData(0)
+											 const float32* _pfProjectionAngles) : configCheckData(0)
 {
 	_clear();
-	_initialize(_iAngleCount, _iDetectorCount, _fDetectorWidth, _pfProjectionAngles,_pfExtraDetectorOffsets);
+	_initialize(_iAngleCount, _iDetectorCount, _fDetectorWidth, _pfProjectionAngles);
 }
 
 //----------------------------------------------------------------------------------------
@@ -156,8 +155,7 @@ bool CProjectionGeometry2D::initialize(const Config& _cfg)
 bool CProjectionGeometry2D::_initialize(int _iProjectionAngleCount, 
 									    int _iDetectorCount, 
 									    float32 _fDetectorWidth, 
-									    const float32* _pfProjectionAngles,
-										const float32* _pfExtraDetectorOffsets)
+									    const float32* _pfProjectionAngles)
 {
 	if (m_bInitialized) {
 		clear();