<commit_before>#include "itkCudaFDKBackProjectionImageFilter.h"

#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIteratorWithIndex.h>
#include <itkLinearInterpolateImageFunction.h>

namespace itk
{

void
CudaFDKBackProjectionImageFilter
::GenerateData()
{
  this->AllocateOutputs();
  this->UpdateAngularWeights();
  
  std::vector<double> angWeights = this->GetAngularWeights();
  
  OutputImageRegionType region = this->GetOutput()->GetRequestedRegion();

  const unsigned int Dimension = ImageType::ImageDimension;
  const unsigned int nProj = this->GetInput(1)->GetLargestPossibleRegion().GetSize(Dimension-1);

  // Ramp factor is the correction for ramp filter which did not account for the divergence of the beam
  const GeometryPointer geometry = dynamic_cast<GeometryType *>(this->GetGeometry().GetPointer());
  double rampFactor = geometry->GetSourceToDetectorDistance() / geometry->GetSourceToIsocenterDistance();
  rampFactor *= 0.5; // Factor 1/2 in eq 176, page 106, Kak & Slaney

  // Rotation center (assumed to be at 0 yet)
  ImageType::PointType rotCenterPoint;
  rotCenterPoint.Fill(0.0);
  ContinuousIndex<double, Dimension> rotCenterIndex;
  this->GetInput(0)->TransformPhysicalPointToContinuousIndex(rotCenterPoint, rotCenterIndex);
  
  // Load dimensions arguments
  int3 vol_dim;
  vol_dim.x = this->GetOutput()->GetLargestPossibleRegion().GetSize()[0];
  vol_dim.y = this->GetOutput()->GetLargestPossibleRegion().GetSize()[1];
  vol_dim.z = this->GetOutput()->GetLargestPossibleRegion().GetSize()[2];

  int2 img_dim;
  img_dim.x = this->GetInput(1)->GetLargestPossibleRegion().GetSize()[0];
  img_dim.y = this->GetInput(1)->GetLargestPossibleRegion().GetSize()[1];

  // Cuda init
  float *dev_vol;
  cudaArray *dev_img;
  float *dev_matrix;
  CUDA_reconstruct_conebeam_init (img_dim, vol_dim, dev_vol, dev_img, dev_matrix);

  // Go over each projection
  for(unsigned int iProj=0; iProj<nProj; iProj++)
    {
    // Extract the current slice
    ProjectionImagePointer projection = this->GetProjection(iProj, angWeights[iProj] * rampFactor);

    // Index to index matrix normalized to have a correct backprojection weight (1 at the isocenter)
    ProjectionMatrixType matrix = GetIndexToIndexProjectionMatrix(iProj, projection);
    double perspFactor = matrix[Dimension-1][Dimension];
    for(unsigned int j=0; j<Dimension; j++)
      perspFactor += matrix[Dimension-1][j] * rotCenterIndex[j];
    matrix /= perspFactor;<commit_msg>Account for index shift + mysterious 0.5 shift in cuda<commit_after>#include "itkCudaFDKBackProjectionImageFilter.h"

#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIteratorWithIndex.h>
#include <itkLinearInterpolateImageFunction.h>

namespace itk
{

void
CudaFDKBackProjectionImageFilter
::GenerateData()
{
  this->AllocateOutputs();
  this->UpdateAngularWeights();
  
  std::vector<double> angWeights = this->GetAngularWeights();
  
  OutputImageRegionType region = this->GetOutput()->GetRequestedRegion();

  const unsigned int Dimension = ImageType::ImageDimension;
  const unsigned int nProj = this->GetInput(1)->GetLargestPossibleRegion().GetSize(Dimension-1);

  // Ramp factor is the correction for ramp filter which did not account for the divergence of the beam
  const GeometryPointer geometry = dynamic_cast<GeometryType *>(this->GetGeometry().GetPointer());
  double rampFactor = geometry->GetSourceToDetectorDistance() / geometry->GetSourceToIsocenterDistance();
  rampFactor *= 0.5; // Factor 1/2 in eq 176, page 106, Kak & Slaney

  // Rotation center (assumed to be at 0 yet)
  ImageType::PointType rotCenterPoint;
  rotCenterPoint.Fill(0.0);
  ContinuousIndex<double, Dimension> rotCenterIndex;
  this->GetInput(0)->TransformPhysicalPointToContinuousIndex(rotCenterPoint, rotCenterIndex);
  
  // Load dimensions arguments
  int3 vol_dim;
  vol_dim.x = this->GetOutput()->GetLargestPossibleRegion().GetSize()[0];
  vol_dim.y = this->GetOutput()->GetLargestPossibleRegion().GetSize()[1];
  vol_dim.z = this->GetOutput()->GetLargestPossibleRegion().GetSize()[2];

  int2 img_dim;
  img_dim.x = this->GetInput(1)->GetLargestPossibleRegion().GetSize()[0];
  img_dim.y = this->GetInput(1)->GetLargestPossibleRegion().GetSize()[1];

  // Cuda init
  float *dev_vol;
  cudaArray *dev_img;
  float *dev_matrix;
  CUDA_reconstruct_conebeam_init (img_dim, vol_dim, dev_vol, dev_img, dev_matrix);

  // Go over each projection
  for(unsigned int iProj=0; iProj<nProj; iProj++)
    {
    // Extract the current slice
    ProjectionImagePointer projection = this->GetProjection(iProj, angWeights[iProj] * rampFactor);

    // Index to index matrix normalized to have a correct backprojection weight (1 at the isocenter)
    ProjectionMatrixType matrix = GetIndexToIndexProjectionMatrix(iProj, projection);

    // We correct the matrix for non zero indexes
    itk::Matrix<double, 4, 4> matrixIdxVol;
    matrixIdxVol.SetIdentity();
    for(unsigned int i=0; i<3; i++)
      matrixIdxVol[i][3] = this->GetOutput()->GetLargestPossibleRegion().GetIndex()[i];
    itk::Matrix<double, 3, 3> matrixIdxProj;
    matrixIdxProj.SetIdentity();
    for(unsigned int i=0; i<2; i++)
      matrixIdxProj[i][2] = -1*(projection->GetLargestPossibleRegion().GetIndex()[i])+0.5; //SR: 0.5 for 2D texture?

    matrix = matrixIdxProj.GetVnlMatrix() * matrix.GetVnlMatrix() * matrixIdxVol.GetVnlMatrix();

    double perspFactor = matrix[Dimension-1][Dimension];
    for(unsigned int j=0; j<Dimension; j++)
      perspFactor += matrix[Dimension-1][j] * rotCenterIndex[j];
    matrix /= perspFactor;
{'fraction_non_alphanumeric': 0.12396990326048012, 'fraction_numerical': 0.013973486205661053, 'mean_word_length': 5.35876993166287, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 66, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3733535', 'n_tokens_mistral': 1847, 'n_tokens_neox': 1784, 'n_words': 380}