ENH: Account for lateral detector offsets in Zeng's forward and back …

…projectors

include/rtkZengBackProjectionImageFilter.hxx

@@ -342,12 +342,14 @@ ZengBackProjectionImageFilter<TInputImage, TOutputImage>::GenerateData()
   indexSlice.Fill(0);
   double dist = NAN, sigmaSlice = NAN;
   double thicknessSlice = this->GetInput(0)->GetSpacing()[2];
-  int    nbProjections = 0;
+  int    nbProjections = indexProj;
   int    startSlice = -1;
   for (auto & angle : list_angle)
   {
     // Get the center of the rotated volume
     m_Transform->SetRotation(0., angle, 0.);
+    m_Transform->SetTranslation(itk::MakeVector(
+      geometry->GetProjectionOffsetsX()[nbProjections], geometry->GetProjectionOffsetsY()[nbProjections], 0.));
     centerRotatedVolume = m_Transform->GetMatrix() * m_centerVolume;
 
     // Set the new origin of the rotate volume according to the center
@@ -357,10 +359,10 @@ ZengBackProjectionImageFilter<TInputImage, TOutputImage>::GenerateData()
     m_ConstantVolumeSource->SetOrigin(originRotatedVolume);
 
     // Set the rotation angle.
-    m_Transform->SetRotation(0., -angle, 0.);
+    m_ResampleImageFilter->SetTransform(m_Transform->GetInverseTransform());
 
     // Extract the projection corresponding to the current angle from the projection stack
-    projRegion.SetIndex(Dimension - 1, nbProjections + indexProj);
+    projRegion.SetIndex(Dimension - 1, nbProjections);
     m_ExtractImageFilter->SetExtractionRegion(projRegion);
     m_ExtractImageFilter->UpdateOutputInformation();
 
@@ -374,12 +376,15 @@ ZengBackProjectionImageFilter<TInputImage, TOutputImage>::GenerateData()
       startSlice += 1;
       indexSlice[Dimension - 1] = startSlice;
       m_ConstantVolumeSource->GetOutput()->TransformIndexToPhysicalPoint(indexSlice, pointSlice);
-      dist = geometry->GetSourceToIsocenterDistances()[nbProjections + indexProj] +
+      dist = geometry->GetSourceToIsocenterDistances()[nbProjections] +
              pointSlice.GetVectorFromOrigin() * m_VectorOrthogonalDetector;
     }
     if (this->GetInput(2))
     {
       m_AttenuationMapTransform->SetRotation(0., angle, 0.);
+      m_Transform->SetTranslation(itk::MakeVector(
+        geometry->GetProjectionOffsetsX()[nbProjections], geometry->GetProjectionOffsetsY()[nbProjections], 0.));
+
       m_AttenuationMapResampleImageFilter->SetOutputOrigin(originRotatedVolume);
       m_AttenuationMapResampleImageFilter->Update();
       rotatedAttenuation = m_AttenuationMapResampleImageFilter->GetOutput();

include/rtkZengForwardProjectionImageFilter.hxx

@@ -327,11 +327,13 @@ ZengForwardProjectionImageFilter<TInputImage, TOutputImage>::GenerateData()
   indexProjection.Fill(0);
   double dist = NAN, sigmaSlice = NAN;
   double thicknessSlice = NAN;
-  int    nbProjections = 0;
+  int    nbProjections = projRegion.GetIndex(Dimension - 1);
   for (auto & angle : list_angle)
   {
     // Set the rotation angle.
     m_Transform->SetRotation(0., angle, 0.);
+    m_Transform->SetTranslation(itk::MakeVector(
+      geometry->GetProjectionOffsetsX()[nbProjections], geometry->GetProjectionOffsetsY()[nbProjections], 0.));
     centerRotatedVolume = m_Transform->GetMatrix() * m_centerVolume;
 
     // Rotate the input volume
@@ -381,7 +383,7 @@ ZengForwardProjectionImageFilter<TInputImage, TOutputImage>::GenerateData()
 
     // Compute the distance between the current slice and the detector
     rotatedVolume->TransformIndexToPhysicalPoint(indexSlice, pointSlice);
-    dist = geometry->GetSourceToIsocenterDistances()[nbProjections + projRegion.GetIndex(Dimension - 1)] +
+    dist = geometry->GetSourceToIsocenterDistances()[nbProjections] +
            pointSlice.GetVectorFromOrigin() * m_VectorOrthogonalDetector;
     int index = 0;
     for (index = nbSlice - 2; index >= 0; index--)
@@ -425,7 +427,7 @@ ZengForwardProjectionImageFilter<TInputImage, TOutputImage>::GenerateData()
     sigmaSlice = pow(m_Alpha * dist + m_SigmaZero, 2.0);
     m_DiscreteGaussianFilter->SetVariance(sigmaSlice);
     // Paste the projection in the output volume
-    indexProjection[Dimension - 1] = nbProjections + projRegion.GetIndex(Dimension - 1);
+    indexProjection[Dimension - 1] = nbProjections;
     m_PasteImageFilter->SetSourceRegion(m_DiscreteGaussianFilter->GetOutput()->GetLargestPossibleRegion());
     m_PasteImageFilter->SetDestinationIndex(indexProjection);
     m_PasteImageFilter->UpdateLargestPossibleRegion();

test/rtkadjointoperatorstest.cxx

@@ -330,7 +330,7 @@ main(int, char **)
   using GeometryType = rtk::ThreeDCircularProjectionGeometry;
   GeometryType::Pointer geometry_parallel = GeometryType::New();
   for (unsigned int noProj = 0; noProj < NumberOfProjectionImages; noProj++)
-    geometry_parallel->AddProjection(500., 0., noProj * 360. / NumberOfProjectionImages);
+    geometry_parallel->AddProjection(500., 0., noProj * 360. / NumberOfProjectionImages, 16., -12.);
 
 
   std::cout << "\n\n******  Zeng Forward projector ******" << std::endl;