fattori_g/reg23Topograms
0
0
Fork 0
Code Issues 1 Pull Requests Actions Packages Projects Releases Wiki Activity

Initial commit.

Browse Source 
Digitally Reconstructed Topogram library and simple test application.
This commit is contained in:
fattori_g
2021-12-17 10:08:47 +01:00
commit c196ae3dbb
5 changed files with 2982 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

40
reg23Topograms/itkDTRrecon/CMakeLists.txt Normal file
View File

@ -0,0 +1,40 @@
cmake_minimum_required(VERSION 2.8.8)
SET(LIB_NAME itkDTRrecon)
find_package(ITK REQUIRED)
include(${ITK_USE_FILE})
INCLUDE_DIRECTORIES(${ITK_INCLUDE_DIRS})
find_package(VTK REQUIRED)
INCLUDE_DIRECTORIES(${VTK_INCLUDE_DIRS})
FIND_PACKAGE(GDCM REQUIRED)
include(${GDCM_USE_FILE})
INCLUDE_DIRECTORIES(${GDCM_INCLUDE_DIRS})
SET(SRCS
itkImageProcessor.cpp
)
SET(HDR
itkImageProcessor.h itkgSiddonJacobsRayCastInterpolateImageFunction.h itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
)
ADD_LIBRARY(${LIB_NAME} ${SRCS} ${HDR})
SET(LINK_LIBS
${LIB_NAME}
${VTK_LIBRARIES}
${ITK_LIBRARIES}
${GDCM_LIBRARIES}
)
TARGET_LINK_LIBRARIES(
${LINK_LIBS}
)
target_include_directories (${LIB_NAME}
PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
)

1933
reg23Topograms/itkDTRrecon/itkImageProcessor.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

273
reg23Topograms/itkDTRrecon/itkImageProcessor.h Normal file
View File

@ -0,0 +1,273 @@
/*
Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
version of incremental ray-tracing algorithm
gfattori 08.11.2021
*/
#ifndef ITKIMAGEPROCESSOR_H
#define ITKIMAGEPROCESSOR_H
#include "itkCommand.h"
#include "itkPowellOptimizer.h"
#include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
#include "itkEuler3DTransform.h"
#include "itkTwoProjectionImageRegistrationMethod.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkFlipImageFilter.h"
#include "itkChangeInformationImageFilter.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h"
#include "itkGDCMImageIO.h"
#include "itkMetaDataObject.h"
#include "gdcmGlobal.h"
#include "itkImageToVTKImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImageSeriesReader.h"
namespace itk
{
class CommandIterationUpdate : public itk::Command
{
public:
using Self = CommandIterationUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
CommandIterationUpdate() = default;
public:
using OptimizerType = itk::PowellOptimizer;
using OptimizerPointer = const OptimizerType *;
void
Execute(itk::Object * caller, const itk::EventObject & event) override
{
Execute((const itk::Object *)caller, event);
}
void
Execute(const itk::Object * object, const itk::EventObject & event) override
{
auto optimizer = dynamic_cast<OptimizerPointer>(object);
if (typeid(event) != typeid(itk::IterationEvent))
{
return;
}
std::cout << "Iteration: " << optimizer->GetCurrentIteration() << std::endl;
std::cout << "Similarity: " << optimizer->GetValue() << std::endl;
std::cout << "Position: " << optimizer->GetCurrentPosition() << std::endl;
}
};
/* reference string required for comparison with tag values */
static const char *ImageOrientationStrings[] = {
"HFS",
"FFS"
};
class ITK_EXPORT itkImageProcessor : public Object
{
constexpr static unsigned int Dimension = 3;
public:
/** Standard "Self" typedef. */
typedef itkImageProcessor Self;
/** Standard "Superclass" typedef. */
typedef Object Superclass;
/** Smart pointer typedef support */
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
/** Method of creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(itkImageProcessor, Object);
/** Input data load methods*/
void load3DVolume(const char *);
int load3DSerie(const char* );
int load2D(const char *);
/** Projection angles - Gantry angle */
void SetProjectionAngle1(double);
void SetProjectionAngle2(double);
/** Source to axis distance - SAD */
void SetSCD(double);
/** Intensity threshold used for image projection,
* any voxel below threshold is ignored */
void SetIntensityThreshold(double);
/** Custom settings of the projection images */
void SetCustom_Isocenter(double,double,double);
void SetCustom_2Dres(double,double,double,double);
void SetCustom_2Dsize(int,int,int,int);
void SetCustom_2Dcenter(double dX1, double dY1, double dX2, double dY2);
/** Set transform parameters for 3D Volume */
void SetInitialTranslations(double, double, double);
void SetInitialRotations(double, double, double);
/** Get the Localizer intensity window and center for rendering */
double* GetImageIntensityWindow(int );
/** 2D-3D registration routines */
// int InitializeRegistration();
void InitializeProjector();
// void RunRegistration();
/** Compute Digital Localizers */
void GetProjectionImages();
/** Conveniency methods to get VTK images for rendering */
vtkImageData* GetProjection1VTK();
vtkImageData* GetProjection2VTK();
vtkImageData* GetLocalizer1VTK();
vtkImageData* GetLocalizer2VTK();
/** Debug writers */
void WriteProjectionImages();
void Write2DImages();
protected:
itkImageProcessor();
virtual ~itkImageProcessor();
void PrintSelf(std::ostream& os, Indent indent) const;
private:
itkImageProcessor(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
void Finalize();
/** Various pixel types */
using InternalPixelType = float;
using PixelType2D = short;
using PixelType3D = short;
using OutputPixelType = unsigned char;
/** Image types */
using ImageType2D = itk::Image<PixelType3D, Dimension>;
using ImageType3D = itk::Image<PixelType3D, Dimension>;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
/** The following lines define each of the components used in the
registration: The transform, optimizer, metric, interpolator and
the registration method itself. */
using InternalImageType = itk::Image<InternalPixelType, Dimension>;
using TransformType = itk::Euler3DTransform<double>;
using OptimizerType = itk::PowellOptimizer;
using MetricType = itk::NormalizedCorrelationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
using RegistrationType = itk::TwoProjectionImageRegistrationMethod<InternalImageType, InternalImageType>;
using ParametersType = RegistrationType::ParametersType;
/** Image reader types */
using ImageReaderType2D = itk::ImageFileReader<ImageType2D>;
using ImageReaderType3D = itk::ImageFileReader<ImageType3D>;
using ImageSeriesReaderType3D = itk::ImageSeriesReader<ImageType3D>;
using ImageIOType = itk::GDCMImageIO;
/** widely used filters: flip, cast, rescale, ... */
using FlipFilterType = itk::FlipImageFilter<InternalImageType>;
using Input2DRescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, InternalImageType>;
using CastFilterType3D = itk::CastImageFilter<ImageType3D, InternalImageType>;
using ChangeInformationFilterType = itk::ChangeInformationImageFilter<InternalImageType>;
using ChangeInformationFilterInputType = itk::ChangeInformationImageFilter<ImageType3D>;
/** ITK to VTK filter */
using ITKtoVTKFilterType = itk::ImageToVTKImageFilter<OutputImageType>;
/** Enum type for Orientation */
enum eImageOrientationType { HFS = 0, FFS = 1};
MetricType::Pointer metric;
TransformType::Pointer transform;
OptimizerType::Pointer optimizer;
InterpolatorType::Pointer interpolator1;
InterpolatorType::Pointer interpolator2;
RegistrationType::Pointer registration;
ImageReaderType2D::Pointer imageReader2D;
ImageReaderType3D::Pointer imageReader3D;
ImageSeriesReaderType3D::Pointer imageSeriesReader3D;
CommandIterationUpdate::Pointer observer;
FlipFilterType::Pointer flipFilter1;
FlipFilterType::Pointer flipFilter2;
FlipFilterType::Pointer flipFilterLocalizer1;
FlipFilterType::Pointer flipFilterLocalizer2;
CastFilterType3D::Pointer caster3D;
CastFilterType3D::Pointer caster2D1;
CastFilterType3D::Pointer caster2D2;
CastFilterType3D::Pointer caster2DInput;
TransformType::InputPointType isocenter3D;
ITKtoVTKFilterType::Pointer toVTK2D1;
ITKtoVTKFilterType::Pointer toVTK2D2;
ITKtoVTKFilterType::Pointer toVTKLocalizer1;
ITKtoVTKFilterType::Pointer toVTKLocalizer2;
ImageType3D::Pointer image3DIn;
InternalImageType::Pointer image2D1In;
InternalImageType::Pointer image2D2In;
Input2DRescaleFilterType::Pointer image2D1Rescaler;
Input2DRescaleFilterType::Pointer image2D2Rescaler;
InternalImageType::Pointer TMPImg1;
InternalImageType::Pointer TMPImg2;
ChangeInformationFilterInputType::Pointer image3DIECFilt;
double image1res[2], image2res[2];
double image1center[2], image2center[2];
int image1Size[2], image2Size[2];
double isocenter[3];
double scd;
double threshold;
double projAngle1, projAngle2;
double TZero[3], RZero[3];
double image3DCOV[3];
bool customized_iso,customized_2DCX,customized_2DRES, customized_2DSIZE;
bool image2D1Loaded, image2D2Loaded, image3DLoaded;
double image1IntensityWindow[2], image2IntensityWindow[2];
};
}
#endif

248
reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h Normal file
View File

@ -0,0 +1,248 @@
/*
Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
version of incremental ray-tracing algorithm
gfattori 08.11.2021
*/
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
/*=========================================================================
Calculate DRR from a CT dataset using incremental ray-tracing algorithm
The algorithm was initially proposed by Robert Siddon and improved by
Filip Jacobs etc.
-------------------------------------------------------------------------
References:
R. L. Siddon, "Fast calculation of the exact radiological path for a
threedimensional CT array," Medical Physics 12, 252-55 (1985).
F. Jacobs, E. Sundermann, B. De Sutter, M. Christiaens, and I. Lemahieu,
"A fast algorithm to calculate the exact radiological path through a pixel
or voxel space," Journal of Computing and Information Technology ?
CIT 6, 89-94 (1998).
=========================================================================*/
#ifndef itkgSiddonJacobsRayCastInterpolateImageFunction_h
#define itkgSiddonJacobsRayCastInterpolateImageFunction_h
#include "itkInterpolateImageFunction.h"
#include "itkTransform.h"
#include "itkVector.h"
#include "itkEuler3DTransform.h"
namespace itk
{
/** \class gSiddonJacobsRayCastInterpolateImageFunction
* \brief Projective interpolation of an image at specified positions.
*
* SiddonJacobsRayCastInterpolateImageFunction casts rays through a 3-dimensional
* image
* \warning This interpolator works for 3-dimensional images only.
*
* \ingroup ImageFunctions
* \ingroup TwoProjectionRegistration
*/
template <typename TInputImage, typename TCoordRep = float>
class gSiddonJacobsRayCastInterpolateImageFunction : public InterpolateImageFunction<TInputImage, TCoordRep>
{
public:
ITK_DISALLOW_COPY_AND_ASSIGN(gSiddonJacobsRayCastInterpolateImageFunction);
/** Standard class type alias. */
using Self = gSiddonJacobsRayCastInterpolateImageFunction;
using Superclass = InterpolateImageFunction<TInputImage, TCoordRep>;
using Pointer = SmartPointer<Self>;
using ConstPointer = SmartPointer<const Self>;
/** Constants for the image dimensions */
static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension;
using TransformType = Euler3DTransform<TCoordRep>;
using TransformPointer = typename TransformType::Pointer;
using InputPointType = typename TransformType::InputPointType;
using OutputPointType = typename TransformType::OutputPointType;
using TransformParametersType = typename TransformType::ParametersType;
using TransformJacobianType = typename TransformType::JacobianType;
using PixelType = typename Superclass::InputPixelType;
using SizeType = typename TInputImage::SizeType;
using DirectionType = Vector<TCoordRep, 3>;
/** Type of the Interpolator Base class */
using InterpolatorType = InterpolateImageFunction<TInputImage, TCoordRep>;
using InterpolatorPointer = typename InterpolatorType::Pointer;
/** Run-time type information (and related methods). */
itkTypeMacro(gSiddonJacobsRayCastInterpolateImageFunction, InterpolateImageFunction);
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** OutputType type alias support. */
using OutputType = typename Superclass::OutputType;
/** InputImageType type alias support. */
using InputImageType = typename Superclass::InputImageType;
/** InputImageConstPointer type alias support. */
using InputImageConstPointer = typename Superclass::InputImageConstPointer;
/** RealType type alias support. */
using RealType = typename Superclass::RealType;
/** Dimension underlying input image. */
static constexpr unsigned int ImageDimension = Superclass::ImageDimension;
/** Point type alias support. */
using PointType = typename Superclass::PointType;
/** Index type alias support. */
using IndexType = typename Superclass::IndexType;
/** ContinuousIndex type alias support. */
using ContinuousIndexType = typename Superclass::ContinuousIndexType;
/** \brief
* Interpolate the image at a point position.
*
* Returns the interpolated image intensity at a
* specified point position. No bounds checking is done.
* The point is assume to lie within the image buffer.
*
* ImageFunction::IsInsideBuffer() can be used to check bounds before
* calling the method.
*/
OutputType
Evaluate(const PointType & point) const override;
/** Interpolate the image at a continuous index position
*
* Returns the interpolated image intensity at a
* specified index position. No bounds checking is done.
* The point is assume to lie within the image buffer.
*
* Subclasses must override this method.
*
* ImageFunction::IsInsideBuffer() can be used to check bounds before
* calling the method.
*/
OutputType
EvaluateAtContinuousIndex(const ContinuousIndexType & index) const override;
virtual void
Initialize();
/** Connect the Transform. */
itkSetObjectMacro(Transform, TransformType);
/** Get a pointer to the Transform. */
itkGetConstObjectMacro(Transform, TransformType);
/** Set and get the focal point to isocenter distance in mm */
itkSetMacro(FocalPointToIsocenterDistance, double);
itkGetMacro(FocalPointToIsocenterDistance, double);
/** Set and get the Lianc grantry rotation angle in radians */
itkSetMacro(ProjectionAngle, double);
itkGetMacro(ProjectionAngle, double);
/** Set and get the Threshold */
itkSetMacro(Threshold, double);
itkGetMacro(Threshold, double);
/** Check if a point is inside the image buffer.
* \warning For efficiency, no validity checking of
* the input image pointer is done. */
inline bool
IsInsideBuffer(const PointType &) const override
{
return true;
}
bool
IsInsideBuffer(const ContinuousIndexType &) const override
{
return true;
}
bool
IsInsideBuffer(const IndexType &) const override
{
return true;
}
#if !defined(ITKV4_COMPATIBILITY)
SizeType
GetRadius() const override
{
const InputImageType * input = this->GetInputImage();
if (!input)
{
itkExceptionMacro("Input image required!");
}
return input->GetLargestPossibleRegion().GetSize();
}
#endif
protected:
gSiddonJacobsRayCastInterpolateImageFunction();
~gSiddonJacobsRayCastInterpolateImageFunction() override = default;
void
PrintSelf(std::ostream & os, Indent indent) const override;
/// Transformation used to calculate the new focal point position
TransformPointer m_Transform; // Displacement of the volume
// Overall inverse transform used to calculate the ray position in the input space
TransformPointer m_InverseTransform;
// The threshold above which voxels along the ray path are integrated
double m_Threshold;
double m_FocalPointToIsocenterDistance; // Focal point to isocenter distance
double m_ProjectionAngle; // Linac gantry rotation angle in radians
private:
void
ComputeInverseTransform() const;
TransformPointer m_GantryRotTransform; // Gantry rotation transform
TransformPointer m_CamShiftTransform; // Camera shift transform camRotTransform
TransformPointer m_CamRotTransform; // Camera rotation transform
TransformPointer m_ComposedTransform; // Composed transform
PointType m_SourcePoint; // Coordinate of the source in the standard Z projection geometry
PointType m_SourceWorld; // Coordinate of the source in the world coordinate system
};
} // namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
# include "itkgSiddonJacobsRayCastInterpolateImageFunction.hxx"
#endif
#endif

488
reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx Normal file
View File

@ -0,0 +1,488 @@
/*
Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
version of incremental ray-tracing algorithm
gfattori 08.11.2021
*/
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
/*=========================================================================
The algorithm was initially proposed by Robert Siddon and improved by
Filip Jacobs etc.
-------------------------------------------------------------------------
References:
R. L. Siddon, "Fast calculation of the exact radiological path for a
threedimensional CT array," Medical Physics 12, 252-55 (1985).
F. Jacobs, E. Sundermann, B. De Sutter, M. Christiaens, and I. Lemahieu,
"A fast algorithm to calculate the exact radiological path through a pixel
or voxel space," Journal of Computing and Information Technology ?
CIT 6, 89-94 (1998).
=========================================================================*/
#ifndef itkgSiddonJacobsRayCastInterpolateImageFunction_hxx
#define itkgSiddonJacobsRayCastInterpolateImageFunction_hxx
#include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
#include "itkMath.h"
#include <cstdlib>
#include <unistd.h>
namespace itk
{
template <typename TInputImage, typename TCoordRep>
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::gSiddonJacobsRayCastInterpolateImageFunction()
{
m_FocalPointToIsocenterDistance = 1000.; // Focal point to isocenter distance in mm.
m_ProjectionAngle = 0.; // Angle in radians betweeen projection central axis and reference axis
m_Threshold = 0.; // Intensity threshold, below which is ignored.
m_SourcePoint[0] = 0.;
m_SourcePoint[1] = 0.;
m_SourcePoint[2] = 0.;
m_InverseTransform = TransformType::New();
m_InverseTransform->SetComputeZYX(true);
m_ComposedTransform = TransformType::New();
m_ComposedTransform->SetComputeZYX(true);
m_GantryRotTransform = TransformType::New();
m_GantryRotTransform->SetComputeZYX(true);
m_GantryRotTransform->SetIdentity();
m_CamShiftTransform = TransformType::New();
m_CamShiftTransform->SetComputeZYX(true);
m_CamShiftTransform->SetIdentity();
m_CamRotTransform = TransformType::New();
m_CamRotTransform->SetComputeZYX(true);
m_CamRotTransform->SetIdentity();
// constant for converting degrees into radians
const float dtr = (atan(1.0) * 4.0) / 180.0;
m_CamRotTransform->SetRotation(dtr * (-90.0), 0.0, 0.0);
m_Threshold = 0;
}
template <typename TInputImage, typename TCoordRep>
void
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::PrintSelf(std::ostream & os, Indent indent) const
{
this->Superclass::PrintSelf(os, indent);
os << indent << "Threshold: " << m_Threshold << std::endl;
os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
}
template <typename TInputImage, typename TCoordRep>
typename gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(const PointType & point) const
{
float rayVector[3];
IndexType cIndex;
PointType drrPixelWorld; // Coordinate of a DRR pixel in the world coordinate system
OutputType pixval;
float firstIntersection[3];
float alphaX1, alphaXN, alphaXmin, alphaXmax;
float alphaY1, alphaYN, alphaYmin, alphaYmax;
float alphaZ1, alphaZN, alphaZmin, alphaZmax;
float alphaMin, alphaMax;
float alphaX, alphaY, alphaZ, alphaCmin, alphaCminPrev;
float alphaUx, alphaUy, alphaUz;
float alphaIntersectionUp[3], alphaIntersectionDown[3];
float d12, value;
float firstIntersectionIndex[3];
int firstIntersectionIndexUp[3], firstIntersectionIndexDown[3];
int iU, jU, kU;
// Min/max values of the output pixel type AND these values
// represented as the output type of the interpolator
const OutputType minOutputValue = itk::NumericTraits<OutputType>::NonpositiveMin();
const OutputType maxOutputValue = itk::NumericTraits<OutputType>::max();
// If the volume was shifted, recalculate the overall inverse transform
unsigned long int interpMTime = this->GetMTime();
unsigned long int vTransformMTime = m_Transform->GetMTime();
if (interpMTime < vTransformMTime)
{
this->ComputeInverseTransform();
// The m_SourceWorld should be computed here to avoid the repeatedly calculation
// for each projection ray. However, we are in a const function, which prohibits
// the modification of class member variables. So the world coordiate of the source
// point is calculated for each ray as below. Performance improvement may be made
// by using a static variable?
// m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
}
//PointType m_SourcePointSliding = m_SourcePoint;
//m_SourcePointSliding[2] = point[2];
PointType SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
// Get ths input pointers
InputImageConstPointer inputPtr = this->GetInputImage();
typename InputImageType::SizeType sizeCT;
typename InputImageType::RegionType regionCT;
typename InputImageType::SpacingType ctPixelSpacing;
typename InputImageType::PointType ctOrigin;
ctPixelSpacing = inputPtr->GetSpacing();
ctOrigin = inputPtr->GetOrigin();
regionCT = inputPtr->GetLargestPossibleRegion();
sizeCT = regionCT.GetSize();
drrPixelWorld = m_InverseTransform->TransformPoint(point);
//slide source along z axis...
SourceWorld[2] = drrPixelWorld[2];
// calculate the detector position for this pixel center by moving
// 2*m_FocalPointToIsocenterDistance from the source in the pixel
// directions
double p1[2], p2[2];
p1[0] = SourceWorld[0];
p1[1] = SourceWorld[1];
p2[0] = drrPixelWorld[0];
p2[1] = drrPixelWorld[1];
double P12D = sqrt( (p2[0]-p1[0]) * (p2[0]-p1[0]) + (p2[1]-p1[1]) * (p2[1]-p1[1]) );
double p12V[2];
p12V[0] = p2[0] - p1[0];
p12V[1] = p2[1] - p1[1];
double p12v [2];
p12v[0] = p12V[0]/P12D;
p12v[1] = p12V[1]/P12D;
double pDest [2];
pDest[0] = p1[0] + m_FocalPointToIsocenterDistance * 2 * p12v[0];
pDest[1] = p1[1] + m_FocalPointToIsocenterDistance * 2 * p12v[1];
// The following is the Siddon-Jacob fast ray-tracing algorithm
//rayVector[0] = drrPixelWorld[0] - SourceWorld[0];
//rayVector[1] = drrPixelWorld[1] - SourceWorld[1];
// correct ray to reach detector
rayVector[0] = pDest[0] - SourceWorld[0];
rayVector[1] = pDest[1] - SourceWorld[1];
rayVector[2] = drrPixelWorld[2] - SourceWorld[2];
/* Calculate the parametric values of the first and the last
intersection points of the ray with the X, Y, and Z-planes that
define the CT volume. */
if (rayVector[0] != 0)
{
alphaX1 = (0.0 - SourceWorld[0]) / rayVector[0];
alphaXN = (sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaXmin = std::min(alphaX1, alphaXN);
alphaXmax = std::max(alphaX1, alphaXN);
}
else
{
alphaXmin = -2;
alphaXmax = 2;
}
if (rayVector[1] != 0)
{
alphaY1 = (0.0 - SourceWorld[1]) / rayVector[1];
alphaYN = (sizeCT[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
alphaYmin = std::min(alphaY1, alphaYN);
alphaYmax = std::max(alphaY1, alphaYN);
}
else
{
alphaYmin = -2;
alphaYmax = 2;
}
if (rayVector[2] != 0)
{
alphaZ1 = (0.0 - SourceWorld[2]) / rayVector[2];
alphaZN = (sizeCT[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaZmin = std::min(alphaZ1, alphaZN);
alphaZmax = std::max(alphaZ1, alphaZN);
}
else
{
alphaZmin = -2;
alphaZmax = 2;
}
/* Get the very first and the last alpha values when the ray
intersects with the CT volume. */
alphaMin = std::max(std::max(alphaXmin, alphaYmin), alphaZmin);
alphaMax = std::min(std::min(alphaXmax, alphaYmax), alphaZmax);
/* Calculate the parametric values of the first intersection point
of the ray with the X, Y, and Z-planes after the ray entered the
CT volume. */
firstIntersection[0] = SourceWorld[0] + alphaMin * rayVector[0];
firstIntersection[1] = SourceWorld[1] + alphaMin * rayVector[1];
firstIntersection[2] = SourceWorld[2] + alphaMin * rayVector[2];
/* Transform world coordinate to the continuous index of the CT volume*/
firstIntersectionIndex[0] = firstIntersection[0] / ctPixelSpacing[0];
firstIntersectionIndex[1] = firstIntersection[1] / ctPixelSpacing[1];
firstIntersectionIndex[2] = firstIntersection[2] / ctPixelSpacing[2];
firstIntersectionIndexUp[0] = (int)ceil(firstIntersectionIndex[0]);
firstIntersectionIndexUp[1] = (int)ceil(firstIntersectionIndex[1]);
firstIntersectionIndexUp[2] = (int)ceil(firstIntersectionIndex[2]);
firstIntersectionIndexDown[0] = (int)floor(firstIntersectionIndex[0]);
firstIntersectionIndexDown[1] = (int)floor(firstIntersectionIndex[1]);
firstIntersectionIndexDown[2] = (int)floor(firstIntersectionIndex[2]);
if (rayVector[0] == 0)
{
alphaX = 2;
}
else
{
alphaIntersectionUp[0] = (firstIntersectionIndexUp[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaIntersectionDown[0] = (firstIntersectionIndexDown[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaX = std::max(alphaIntersectionUp[0], alphaIntersectionDown[0]);
}
if (rayVector[1] == 0)
{
alphaY = 2;
}
else
{
alphaIntersectionUp[1] = (firstIntersectionIndexUp[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
alphaIntersectionDown[1] = (firstIntersectionIndexDown[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
alphaY = std::max(alphaIntersectionUp[1], alphaIntersectionDown[1]);
}
if (rayVector[2] == 0)
{
alphaZ = 2;
}
else
{
alphaIntersectionUp[2] = (firstIntersectionIndexUp[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaIntersectionDown[2] = (firstIntersectionIndexDown[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaZ = std::max(alphaIntersectionUp[2], alphaIntersectionDown[2]);
}
/* Calculate alpha incremental values when the ray intercepts with x, y, and z-planes */
if (rayVector[0] != 0)
{
alphaUx = ctPixelSpacing[0] / fabs(rayVector[0]);
}
else
{
alphaUx = 999;
}
if (rayVector[1] != 0)
{
alphaUy = ctPixelSpacing[1] / fabs(rayVector[1]);
}
else
{
alphaUy = 999;
}
if (rayVector[2] != 0)
{
alphaUz = ctPixelSpacing[2] / fabs(rayVector[2]);
}
else
{
alphaUz = 999;
}
/* Calculate voxel index incremental values along the ray path. */
if (SourceWorld[0] < drrPixelWorld[0])
{
iU = 1;
}
else
{
iU = -1;
}
if (SourceWorld[1] < drrPixelWorld[1])
{
jU = 1;
}
else
{
jU = -1;
}
if (SourceWorld[2] < drrPixelWorld[2])
{
kU = 1;
}
else
{
kU = -1;
}
d12 = 0.0; /* Initialize the sum of the voxel intensities along the ray path to zero. */
/* Initialize the current ray position. */
alphaCmin = std::min(std::min(alphaX, alphaY), alphaZ);
/* Initialize the current voxel index. */
cIndex[0] = firstIntersectionIndexDown[0];
cIndex[1] = firstIntersectionIndexDown[1];
cIndex[2] = firstIntersectionIndexDown[2];
while (alphaCmin < alphaMax) /* Check if the ray is still in the CT volume */
{
/* Store the current ray position */
alphaCminPrev = alphaCmin;
if ((alphaX <= alphaY) && (alphaX <= alphaZ))
{
/* Current ray front intercepts with x-plane. Update alphaX. */
alphaCmin = alphaX;
cIndex[0] = cIndex[0] + iU;
alphaX = alphaX + alphaUx;
}
else if ((alphaY <= alphaX) && (alphaY <= alphaZ))
{
/* Current ray front intercepts with y-plane. Update alphaY. */
alphaCmin = alphaY;
cIndex[1] = cIndex[1] + jU;
alphaY = alphaY + alphaUy;
}
else
{
/* Current ray front intercepts with z-plane. Update alphaZ. */
alphaCmin = alphaZ;
cIndex[2] = cIndex[2] + kU;
alphaZ = alphaZ + alphaUz;
}
if ((cIndex[0] >= 0) && (cIndex[0] < static_cast<IndexValueType>(sizeCT[0])) && (cIndex[1] >= 0) &&
(cIndex[1] < static_cast<IndexValueType>(sizeCT[1])) && (cIndex[2] >= 0) &&
(cIndex[2] < static_cast<IndexValueType>(sizeCT[2])))
{
/* If it is a valid index, get the voxel intensity. */
value = static_cast<float>(inputPtr->GetPixel(cIndex));
if (value > m_Threshold) /* Ignore voxels whose intensities are below the threshold. */
{
d12 += (alphaCmin - alphaCminPrev) * (value - m_Threshold);
}
}
}
if (d12 < minOutputValue)
{
pixval = minOutputValue;
}
else if (d12 > maxOutputValue)
{
pixval = maxOutputValue;
}
else
{
pixval = static_cast<OutputType>(d12);
}
return (pixval);
}
template <typename TInputImage, typename TCoordRep>
typename gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::EvaluateAtContinuousIndex(
const ContinuousIndexType & index) const
{
OutputPointType point;
this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
return this->Evaluate(point);
}
template <typename TInputImage, typename TCoordRep>
void
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::ComputeInverseTransform() const
{
m_ComposedTransform->SetIdentity();
m_ComposedTransform->Compose(m_Transform, 0);
typename TransformType::InputPointType isocenter;
isocenter = m_Transform->GetCenter();
// An Euler 3D transform is used to rotate the volume to simulate the roation of the linac gantry.
// The rotation is about z-axis. After the transform, a AP projection geometry (projecting
// towards positive y direction) is established.
m_GantryRotTransform->SetRotation(0.0, 0.0, -m_ProjectionAngle);
m_GantryRotTransform->SetCenter(isocenter);
m_ComposedTransform->Compose(m_GantryRotTransform, 0);
// An Euler 3D transfrom is used to shift the source to the origin.
typename TransformType::OutputVectorType focalpointtranslation;
focalpointtranslation[0] = -isocenter[0];
focalpointtranslation[1] = m_FocalPointToIsocenterDistance - isocenter[1];
focalpointtranslation[2] = -isocenter[2];
m_CamShiftTransform->SetTranslation(focalpointtranslation);
m_ComposedTransform->Compose(m_CamShiftTransform, 0);
// A Euler 3D transform is used to establish the standard negative z-axis projection geometry. (By
// default, the camera is situated at the origin, points down the negative z-axis, and has an up-
// vector of (0, 1, 0).)
m_ComposedTransform->Compose(m_CamRotTransform, 0);
// The overall inverse transform is computed. The inverse transform will be used by the interpolation
// procedure.
m_ComposedTransform->GetInverse(m_InverseTransform);
this->Modified();
}
template <typename TInputImage, typename TCoordRep>
void
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Initialize()
{
this->ComputeInverseTransform();
m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
}
} // namespace itk
#endif