C++ ASTRA_ASSERT函数代码示例

//----------------------------------------------------------------------------------------
// Iterate
void CCudaFDKAlgorithm3D::run(int _iNrIterations)
{
	// check initialized
	ASTRA_ASSERT(m_bIsInitialized);

	CFloat32ProjectionData3D* pSinoMem = dynamic_cast<CFloat32ProjectionData3D*>(m_pSinogram);
	ASTRA_ASSERT(pSinoMem);
	CFloat32VolumeData3D* pReconMem = dynamic_cast<CFloat32VolumeData3D*>(m_pReconstruction);
	ASTRA_ASSERT(pReconMem);

	const float *filter = NULL;
	if (m_iFilterDataId != -1) {
		const CFloat32ProjectionData2D *pFilterData = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(m_iFilterDataId));
		if (pFilterData)
			filter = pFilterData->getDataConst();
	}

#if 0
	bool ok = true;
	
	ok = astraCudaFDK(pReconMem->getData(), pSinoMem->getDataConst(),
	                  &volgeom, conegeom,
	                  m_bShortScan, m_iGPUIndex, m_iVoxelSuperSampling, filter);

	ASTRA_ASSERT(ok);
#endif

	CCompositeGeometryManager cgm;

	cgm.doFDK(m_pProjector, pReconMem, pSinoMem, m_bShortScan, filter);



}

void CConeProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ,
                                                 int iAngleIndex,
                                                 float32 &fU, float32 &fV) const
{
	ASTRA_ASSERT(iAngleIndex >= 0);
	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);

	float alpha = m_pfProjectionAngles[iAngleIndex];

	// Project point onto optical axis

	// Projector direction is (cos(alpha), sin(alpha))
	// 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;

	// 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;

	// Scale fS to detector plane
	fU = detectorOffsetXToColIndexFloat( (fS * (m_fOriginSourceDistance + m_fOriginDetectorDistance)) / fD );

	ASTRA_DEBUG("alpha: %f, D: %f, V: %f, S: %f, U: %f", alpha, fD, fV, fS, fU);

}

//----------------------------------------------------------------------------------------
// Initializes an instance of the CFloat32Data2D class, allocating (but not initializing) the data block.
bool CFloat32Data3DMemory::_initialize(int _iWidth, int _iHeight, int _iDepth)
{
	// basic checks
	ASTRA_ASSERT(_iWidth > 0);
	ASTRA_ASSERT(_iHeight > 0);
	ASTRA_ASSERT(_iDepth > 0);

	if (m_bInitialized)
	{
		_unInit();
	}
	
	// calculate size
	m_iWidth = _iWidth;
	m_iHeight = _iHeight;
	m_iDepth = _iDepth;
	m_iSize = (size_t)m_iWidth * m_iHeight * m_iDepth;

	// allocate memory for the data, but do not fill it
	m_pfData = NULL;
	m_pCustomMemory = 0;
	_allocateData();

	// initialization complete
	return true;

}

//----------------------------------------------------------------------------------------
// Iterate
void CCudaFDKAlgorithm3D::run(int _iNrIterations)
{
	// check initialized
	ASTRA_ASSERT(m_bIsInitialized);

	const CProjectionGeometry3D* projgeom = m_pSinogram->getGeometry();
	const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(projgeom);
	const CVolumeGeometry3D& volgeom = *m_pReconstruction->getGeometry();

	ASTRA_ASSERT(conegeom);

	CFloat32ProjectionData3DMemory* pSinoMem = dynamic_cast<CFloat32ProjectionData3DMemory*>(m_pSinogram);
	ASTRA_ASSERT(pSinoMem);
	CFloat32VolumeData3DMemory* pReconMem = dynamic_cast<CFloat32VolumeData3DMemory*>(m_pReconstruction);
	ASTRA_ASSERT(pReconMem);


	bool ok = true;

	ok = astraCudaFDK(pReconMem->getData(), pSinoMem->getDataConst(),
	                  &volgeom, conegeom,
	                  m_bShortScan, m_iGPUIndex, m_iVoxelSuperSampling);

	ASTRA_ASSERT(ok);

}

//----------------------------------------------------------------------------------------
// Iterate
void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations)
{
	// check initialized
	ASTRA_ASSERT(m_bIsInitialized);

	CFloat32ProjectionData3DMemory* pSinoMem = dynamic_cast<CFloat32ProjectionData3DMemory*>(m_pSinogram);
	ASTRA_ASSERT(pSinoMem);
	CFloat32VolumeData3DMemory* pReconMem = dynamic_cast<CFloat32VolumeData3DMemory*>(m_pReconstruction);
	ASTRA_ASSERT(pReconMem);

	const CProjectionGeometry3D* projgeom = pSinoMem->getGeometry();
	const CVolumeGeometry3D& volgeom = *pReconMem->getGeometry();

	if (m_bSIRTWeighting) {
		astraCudaBP_SIRTWeighted(pReconMem->getData(),
		                         pSinoMem->getDataConst(),
		                         &volgeom, projgeom,
		                         m_iGPUIndex, m_iVoxelSuperSampling);
	} else {

#if 1
		CCompositeGeometryManager cgm;

		cgm.doBP(m_pProjector, pReconMem, pSinoMem);
#else
		astraCudaBP(pReconMem->getData(), pSinoMem->getDataConst(),
		            &volgeom, projgeom,
		            m_iGPUIndex, m_iVoxelSuperSampling);
#endif
	}

}

//---------------------------------------------------------------------------------------
// Initialize - Config
bool CSirtAlgorithm::initialize(const Config& _cfg)
{
	ASTRA_ASSERT(_cfg.self);
	ConfigStackCheck<CAlgorithm> CC("SirtAlgorithm", this, _cfg);

	// if already initialized, clear first
	if (m_bIsInitialized) {
		clear();
	}

	// initialization of parent class
	if (!CSartAlgorithm::initialize(_cfg)) {
		return false;
	}

	//// init data objects and data projectors
	//_init();

	//// Alpha
	//m_fAlpha = _cfg.self.getOptionNumerical("Alpha", m_fAlpha);
	//CC.markOptionParsed("Alpha");

	// success
	m_bIsInitialized = _check();
	return m_bIsInitialized;
}

//---------------------------------------------------------------------------------------
// Initialize - Config
bool CConeVecProjectionGeometry3D::initialize(const Config& _cfg)
{
	ASTRA_ASSERT(_cfg.self);
	ConfigStackCheck<CProjectionGeometry3D> CC("ConeVecProjectionGeometry3D", this, _cfg);	

	XMLNode* node;

	// TODO: Fix up class hierarchy... this class doesn't fit very well.
	// initialization of parent class
	//CProjectionGeometry3D::initialize(_cfg);

	// Required: DetectorRowCount
	node = _cfg.self->getSingleNode("DetectorRowCount");
	ASTRA_CONFIG_CHECK(node, "ConeVecProjectionGeometry3D", "No DetectorRowCount tag specified.");
	m_iDetectorRowCount = boost::lexical_cast<int>(node->getContent());
	ASTRA_DELETE(node);
	CC.markNodeParsed("DetectorRowCount");

	// Required: DetectorColCount
	node = _cfg.self->getSingleNode("DetectorColCount");
	ASTRA_CONFIG_CHECK(node, "ConeVecProjectionGeometry3D", "No DetectorColCount tag specified.");
	m_iDetectorColCount = boost::lexical_cast<int>(node->getContent());
	m_iDetectorTotCount = m_iDetectorRowCount * m_iDetectorColCount;
	ASTRA_DELETE(node);
	CC.markNodeParsed("DetectorColCount");

	// Required: Vectors
	node = _cfg.self->getSingleNode("Vectors");
	ASTRA_CONFIG_CHECK(node, "ConeVecProjectionGeometry3D", "No Vectors tag specified.");
	vector<double> data = node->getContentNumericalArrayDouble();
	CC.markNodeParsed("Vectors");
	ASTRA_DELETE(node);
	ASTRA_CONFIG_CHECK(data.size() % 12 == 0, "ConeVecProjectionGeometry3D", "Vectors doesn't consist of 12-tuples.");
	m_iProjectionAngleCount = data.size() / 12;
	m_pProjectionAngles = new SConeProjection[m_iProjectionAngleCount];

	for (int i = 0; i < m_iProjectionAngleCount; ++i) {
		SConeProjection& p = m_pProjectionAngles[i];
		p.fSrcX  = data[12*i +  0];
		p.fSrcY  = data[12*i +  1];
		p.fSrcZ  = data[12*i +  2];
		p.fDetUX = data[12*i +  6];
		p.fDetUY = data[12*i +  7];
		p.fDetUZ = data[12*i +  8];
		p.fDetVX = data[12*i +  9];
		p.fDetVY = data[12*i + 10];
		p.fDetVZ = data[12*i + 11];

		// The backend code currently expects the corner of the detector, while
		// the matlab interface supplies the center
		p.fDetSX = data[12*i +  3] - 0.5f * m_iDetectorRowCount * p.fDetVX - 0.5f * m_iDetectorColCount * p.fDetUX;
		p.fDetSY = data[12*i +  4] - 0.5f * m_iDetectorRowCount * p.fDetVY - 0.5f * m_iDetectorColCount * p.fDetUY;
		p.fDetSZ = data[12*i +  5] - 0.5f * m_iDetectorRowCount * p.fDetVZ - 0.5f * m_iDetectorColCount * p.fDetUZ;
	}

	// success
	m_bInitialized = _check();
	return m_bInitialized;
}

//---------------------------------------------------------------------------------------
// Initialize - Config
bool CCudaBackProjectionAlgorithm3D::initialize(const Config& _cfg)
{
	ASTRA_ASSERT(_cfg.self);
	ConfigStackCheck<CAlgorithm> CC("CudaBackProjectionAlgorithm3D", this, _cfg);	

	// if already initialized, clear first
	if (m_bIsInitialized) {
		clear();
	}

	// initialization of parent class
	if (!CReconstructionAlgorithm3D::initialize(_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");


	m_iVoxelSuperSampling = 1;
	if (pCudaProjector)
		m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling();
	m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling);
	CC.markOptionParsed("VoxelSuperSampling");

	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");
	}

	// success
	m_bIsInitialized = _check();
	return m_bIsInitialized;
}

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;
}

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);
}

//----------------------------------------------------------------------------------------
// Initialization with a Config object
bool CProjectionGeometry2D::initialize(const Config& _cfg)
{
	ASTRA_ASSERT(_cfg.self);
	ConfigStackCheck<CProjectionGeometry2D> CC("ProjectionGeometry2D", this, _cfg);	

	// uninitialize if the object was initialized before
	if (m_bInitialized)	{
		clear();
	}

	// Required: DetectorWidth
	XMLNode* node = _cfg.self->getSingleNode("DetectorWidth");
	ASTRA_CONFIG_CHECK(node, "ProjectionGeometry2D", "No DetectorWidth tag specified.");
	m_fDetectorWidth = boost::lexical_cast<float32>(node->getContent());
	ASTRA_DELETE(node);
	CC.markNodeParsed("DetectorWidth");

	// Required: DetectorCount
	node = _cfg.self->getSingleNode("DetectorCount");
	ASTRA_CONFIG_CHECK(node, "ProjectionGeometry2D", "No DetectorCount tag specified.");
	m_iDetectorCount = boost::lexical_cast<int>(node->getContent());
	ASTRA_DELETE(node);
	CC.markNodeParsed("DetectorCount");

	// Required: ProjectionAngles
	node = _cfg.self->getSingleNode("ProjectionAngles");
	ASTRA_CONFIG_CHECK(node, "ProjectionGeometry2D", "No ProjectionAngles tag specified.");
	vector<float32> angles = node->getContentNumericalArray();
	delete node;
	m_iProjectionAngleCount = angles.size();
	ASTRA_CONFIG_CHECK(m_iProjectionAngleCount > 0, "ProjectionGeometry2D", "Not enough ProjectionAngles specified.");
	m_pfProjectionAngles = new float32[m_iProjectionAngleCount];
	for (int i = 0; i < m_iProjectionAngleCount; i++) {
		m_pfProjectionAngles[i] = angles[i];
	}
	CC.markNodeParsed("ProjectionAngles");

	vector<float32> offset = _cfg.self->getOptionNumericalArray("ExtraDetectorOffset");
	m_pfExtraDetectorOffset = new float32[m_iProjectionAngleCount];
	if (offset.size() == (size_t)m_iProjectionAngleCount) {
		for (int i = 0; i < m_iProjectionAngleCount; i++) {
			m_pfExtraDetectorOffset[i] = offset[i];
		}
	} else {
		for (int i = 0; i < m_iProjectionAngleCount; i++) {
			m_pfExtraDetectorOffset[i] = 0.0f;
		}	
	}
	CC.markOptionParsed("ExtraDetectorOffset");

	// some checks
	ASTRA_CONFIG_CHECK(m_iDetectorCount > 0, "ProjectionGeometry2D", "DetectorCount should be positive.");
	ASTRA_CONFIG_CHECK(m_fDetectorWidth > 0.0f, "ProjectionGeometry2D", "DetectorWidth should be positive.");
	ASTRA_CONFIG_CHECK(m_pfProjectionAngles != NULL, "ProjectionGeometry2D", "ProjectionAngles not initialized");

	// Interface class, so don't return true
	return false;
}

CVector3D CFanFlatVecProjectionGeometry2D::getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex /* = 0 */)
{
	CVector3D vOutput(0.0f, 0.0f, 0.0f);

	// not implemented
	ASTRA_ASSERT(false);

	return vOutput;
}

//----------------------------------------------------------------------------------------
CVector3D CSparseMatrixProjectionGeometry2D::getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex)
{
	CVector3D vOutput(0.0f, 0.0f, 0.0f);

	// not implemented, yet
	ASTRA_ASSERT(false);

	return vOutput;
}

void CConeVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ,
                                                 int iAngleIndex,
                                                 float32 &fU, float32 &fV) const
{
	ASTRA_ASSERT(iAngleIndex >= 0);
	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);

	double fUX, fUY, fUZ, fUC;
	double fVX, fVY, fVZ, fVC;
	double fDX, fDY, fDZ, fDC;

	computeBP_UV_Coeffs(m_pProjectionAngles[iAngleIndex],
	                    fUX, fUY, fUZ, fUC, fVX, fVY, fVZ, fVC, fDX, fDY, fDZ, fDC);

	// The -0.5f shifts from corner to center of detector pixels
	double fD = fDX*fX + fDY*fY + fDZ*fZ + fDC;
	fU = (fUX*fX + fUY*fY + fUZ*fZ + fUC) / fD - 0.5f;
	fV = (fVX*fX + fVY*fY + fVZ*fZ + fVC) / fD - 0.5f;
}

//---------------------------------------------------------------------------------------
// Initialize - Config
bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg)
{
	ASTRA_ASSERT(_cfg.self);
	ConfigStackCheck<CAlgorithm> CC("CudaForwardProjectionAlgorithm", this, _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");
		}
	}
	CC.markNodeParsed("ProjectorId");


	
	// sinogram data
	node = _cfg.self.getSingleNode("ProjectionDataId");
	ASTRA_CONFIG_CHECK(node, "FP_CUDA", "No ProjectionDataId tag specified.");
	int id = boost::lexical_cast<int>(node.getContent());
	m_pSinogram = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id));
	CC.markNodeParsed("ProjectionDataId");

	// volume data
	node = _cfg.self.getSingleNode("VolumeDataId");
	ASTRA_CONFIG_CHECK(node, "FP_CUDA", "No VolumeDataId tag specified.");
	id = boost::lexical_cast<int>(node.getContent());
	m_pVolume = dynamic_cast<CFloat32VolumeData2D*>(CData2DManager::getSingleton().get(id));
	CC.markNodeParsed("VolumeDataId");

	// 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");

	// 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
	return check();
}