#include "clitkBackProjectImageFilter.h"
#include "itkContinuousIndex.h"
#include "vnl/vnl_math.h"
+#include "itkLinearInterpolateImageFunction.h"
namespace clitk
{
this->m_SourceToAxis = 1000.0;
this->m_EdgePaddingValue = itk::NumericTraits<OutputPixelType>::Zero;//density images
this->m_RigidTransformMatrix.SetIdentity();
+ this->m_PanelShift[0] = 0.;
+ this->m_PanelShift[1] = 0.;
//Parameters for output
this->m_OutputSpacing.Fill(1);
template <class InputImageType, class OutputImageType>
void
BackProjectImageFilter<InputImageType, OutputImageType>
+#if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR > 12 ) || ( ITK_VERSION_MAJOR > 4 ))
+ ::Initialize( void )
+#else
::Initialize( void ) throw (itk::ExceptionObject)
+#endif
{
//Change the origin of the 2D input
typename InputImageType::ConstPointer inputPtr=this->GetInput();
BackProjectImageFilter<InputImageType, OutputImageType>
::CalculateProjectionMatrix( void )
{
- InputSpacingType inputSpacing=this->GetInput()->GetSpacing();
+ //InputSpacingType inputSpacing=this->GetInput()->GetSpacing();
// Projection on YZ plane+pixelTrans
itk::Matrix<double,3,4> temp;
//-----------------------------------------------------------------------
template <class InputImageType, class OutputImageType>
void
- BackProjectImageFilter<InputImageType, OutputImageType>
- ::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread, int threadId )
+ BackProjectImageFilter<InputImageType, OutputImageType>::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread, itk::ThreadIdType threadId )
{
//Projection pointer
InputImageConstPointer inputPtr=this->GetInput();
- InputPixelType * beginPtr=const_cast<InputPixelType *>(this->GetInput()->GetBufferPointer());
- InputPixelType * pp;
//Volume pointer
OutputImagePointer outputPTr= this->GetOutput();
iPoint.Fill(itk::NumericTraits<double>::Zero);
OutputIndexType oIndex;
ContinuousInputIndexType iIndex;
- InputSizeType inputSize=inputPtr->GetLargestPossibleRegion().GetSize();
- double dx,dy,dxm,dym;
- int lx, ly;
//Get the first output coordinate
oIndex=iterator.GetIndex();//costly but only once a thread
//Compute the first input coordinate (invert Y/X)
homInputPoint= (m_ProjectionMatrix * homOutputPoint);
- iPoint[0]=-homInputPoint[0]/homInputPoint[2];
- iPoint[1]=homInputPoint[1]/homInputPoint[2];
-
+ iPoint[0]=-homInputPoint[0]/homInputPoint[2] + m_PanelShift[0];
+ iPoint[1]=homInputPoint[1]/homInputPoint[2] + m_PanelShift[1];
+
+ typedef itk::LinearInterpolateImageFunction< InputImageType, double > InterpolatorType;
+ typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
+ interpolator->SetInputImage(this->GetInput());
+
//Run over all output voxels
for (unsigned int i=0; i<outputSizeForThread[2]; i++)
{
{
for (unsigned int k=0; k<outputSizeForThread[0]; k++)
{
- iPoint[0]=homInputPoint[0]/homInputPoint[2];
- iPoint[1]=homInputPoint[1]/homInputPoint[2];
+ iPoint[0]=-homInputPoint[0]/homInputPoint[2] + m_PanelShift[0];
+ iPoint[1]=homInputPoint[1]/homInputPoint[2] + m_PanelShift[1];
//Check wether inside, convert to index (use modified with correct origin)
- if( m_ModifiedInput->TransformPhysicalPointToContinuousIndex(iPoint, iIndex) )
- {
- //Own (fast bilinear) interpolation
- lx = (int)floor(iIndex[0]); dx = iIndex[0]-lx; dxm = 1.-dx;
- ly = (int)floor(iIndex[1]); dy = iIndex[1]-ly; dym = 1.-dy;
- pp = beginPtr + ly*inputSize[0]+lx;
- value =static_cast<OutputPixelType>( ( dxm * dym*(double)(*pp)
- + dx * dym*(double)(*(pp+1))
- + dxm* dy *(double)(*(pp + inputSize[0]))
- + dx * dy *(double)(*(pp + inputSize[0]+1))) );
-
- }
+ if (m_ModifiedInput->TransformPhysicalPointToContinuousIndex(iPoint, iIndex) && interpolator->IsInsideBuffer(iIndex))
+ value = interpolator->EvaluateAtContinuousIndex(iIndex);
//Outside: padding value
- else value=m_EdgePaddingValue;
-
+ else
+ value=m_EdgePaddingValue;
//Set it
iterator.Set(value);