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

Merge branch 'Siddon' into 'ScoutRTUIDevel'

Browse Source 
Siddon

See merge request cpt_bioeng/drt!31
This commit is contained in:
fattori_g
2023-06-06 15:21:24 +02:00
parent 1c247296b2 39a9c1e47b
commit a59de4370c
6 changed files with 499 additions and 421 deletions
Download Patch File Download Diff File Expand all files Collapse all files

100
reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
View File

@ -70,17 +70,11 @@ itkImageProcessor::itkImageProcessor()
}
}
m_3DInputChangeInformationToZero =
ChangeInformationFilterType::New();
m_Projection1VTKTransform = vtkMatrix4x4::New();
m_Projection1VTKTransform->Identity();
m_Projection2VTKTransform = vtkMatrix4x4::New();
m_Projection2VTKTransform->Identity();
/* Set to NULL the metainfo */
m_CTMetaInfo = NULL;
@ -230,6 +224,38 @@ double itkImageProcessor::GetSCD2(){
m_DRTImage2MetaInfo ->GetSCD();
}
double itkImageProcessor::GetPanelOffsetPGeo(tProjOrientationType ImgPrj){
if(m_CTMetaInfo == NULL ||
m_DRTGeometryMetaInfo == NULL)
return 0.;
switch (ImgPrj) {
case tProjOrientationType::PA:
return
this->CalcPanelOffsetPGeo(m_CTMetaInfo->GetPatientOrientation(),
-m_DRTGeometryMetaInfo->GetPanel1Offset());
break;
case tProjOrientationType::LAT:
return
this->CalcPanelOffsetPGeo(m_CTMetaInfo->GetPatientOrientation(),
-m_DRTGeometryMetaInfo->GetPanel2Offset());
break;
case tProjOrientationType::NA:
return 0.;
break;
default:
return 0.;
break;
}
}
void itkImageProcessor::SetPanelOffsets(double dOff1, double dOff2)
{
m_DRTGeometryMetaInfo->SetPanel1Offset(dOff1);
@ -551,6 +577,7 @@ int itkImageProcessor::fill3DVolumeMeta(
m_CTMetaInfo->SetSpacing(m_InputImage->GetSpacing());
m_CTMetaInfo->SetSize(
m_InputImage->GetBufferedRegion().GetSize() );
std::cout<<"ORIGIN: "<<m_InputImage->GetOrigin()<<std::endl;
m_CTMetaInfo->SetOriginLPS(m_InputImage->GetOrigin());
m_CTMetaInfo->SetImageDirections(m_InputImage->GetDirection());
ImageType3D::PointType
@ -612,25 +639,8 @@ int itkImageProcessor::fill3DVolumeMeta(
// to be (0,0,0). Because we align the CT isocenter with the central axis, the projection
// geometry is fully defined. The origin of the CT image becomes irrelavent.
ImageType3D::PointType pZeroOrigin;
pZeroOrigin[0] = 0.;
pZeroOrigin[1] = 0.;
pZeroOrigin[2] = 0.;
// this is done internally in the dSiddonFilter.
if (m_3DInputChangeInformationToZero) {
m_3DInputChangeInformationToZero = NULL;
m_3DInputChangeInformationToZero = ChangeInformationFilterType::New();
}
m_3DInputChangeInformationToZero->SetInput(
m_VolumeSourceDupli->GetOutput());
m_3DInputChangeInformationToZero->SetOutputOrigin(
pZeroOrigin);
m_3DInputChangeInformationToZero->UpdateOutputInformation();
m_3DInputChangeInformationToZero->Update();
image3DIn =
m_3DInputChangeInformationToZero->GetOutput();
return EXIT_SUCCESS;
@ -923,7 +933,7 @@ int itkImageProcessor::load2D(const char * pcFName){
ImageDirectionType3D LocalizerImagDirectionDCM =
imageReader2D->GetOutput()->GetDirection();
// std::cout<<"LocalizerImagDirectionDCM " <<LocalizerImagDirectionDCM <<std::endl;
// std::cout<<"LocalizerImagDirectionDCM " <<LocalizerImagDirectionDCM <<std::endl;
InternalImageType::DirectionType toIECDirection;
@ -942,11 +952,17 @@ int itkImageProcessor::load2D(const char * pcFName){
}
}
// std::cout<<"toIECDirection " <<toIECDirection <<std::endl;
/* calculate image orientation with respect to fixed IEC */
InternalImageType::DirectionType LocalizerImagDirectionIEC =
toIECDirection * LocalizerImagDirectionDCM;
// std::cout<<"toIECDirection " <<toIECDirection <<std::endl;
// std::cout<<"LocalizerImagDirectionDCM " <<LocalizerImagDirectionDCM <<std::endl;
// std::cout<<"LocalizerImagDirectionIEC " <<LocalizerImagDirectionIEC <<std::endl;
/* we calculate dot products between the Z components */
double dSumPA = 0;
dSumPA =
@ -960,6 +976,8 @@ int itkImageProcessor::load2D(const char * pcFName){
(LocalizerImagDirectionIEC[1][2] * LATElementsIEC[5]) +
(LocalizerImagDirectionIEC[2][2] * LATElementsIEC[8]);
// std::cout<<"dSumPA " <<dSumPA <<std::endl;
// std::cout<<"dSumLAT " <<dSumLAT <<std::endl;
tProjOrientationType
currProjOrientation = NA;
@ -1324,7 +1342,6 @@ void itkImageProcessor::InitializeProjector()
transform1->SetCenter(
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() );
// 2D Image 1
interpolator1->SetProjectionAngle(
dtr * m_DRTImage1MetaInfo->GetProjectionAngleLPS() );
@ -1357,10 +1374,8 @@ void itkImageProcessor::InitializeProjector()
CurrTransform->GetAngleZ()
);
transform2->SetCenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() );
//transform2->Print(std::cout);
// 2D Image 2
interpolator2->SetProjectionAngle(
@ -1378,7 +1393,11 @@ void itkImageProcessor::InitializeProjector()
resampleFilter1 = ResampleFilterType::New();
resampleFilter1->SetInput(
image3DIn);
m_VolumeSourceDupli->GetOutput());
ImageType3D::RegionType lagerReg =
m_VolumeSourceDupli->GetOutput()->GetLargestPossibleRegion();
resampleFilter1->SetDefaultPixelValue(0);
resampleFilter1->SetNumberOfWorkUnits(iNWUnits);
@ -1397,7 +1416,8 @@ void itkImageProcessor::InitializeProjector()
resampleFilter2 = ResampleFilterType::New();
resampleFilter2->SetInput(
image3DIn);
m_VolumeSourceDupli->GetOutput());
resampleFilter2->SetDefaultPixelValue(0);
resampleFilter2->SetNumberOfWorkUnits(iNWUnits);
@ -1595,10 +1615,10 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
m_DRTGeometryMetaInfo->GetProjectionCenterOffset1());
// std::cout<< "///////////////// PGEOM META BEG ///////////////" <<std::endl;
// std::cout<<"NominalIsocenter IEC"<< m_DRTGeometryMetaInfo->GetProjectionCenter() <<std::endl;
// std::cout<<"CALIBRATION NominalIsocenter_wZcorrectionLPS "<<NominalIsocenter_wZcorrectionLPS<<std::endl;
// std::cout<<"CALIBRATION NominalIsocenterZero_wZcorrectionLPS "<<NominalIsocenterZero_wZcorrectionLPS<<std::endl;
// std::cout<<"CALIBRATION IsocenterOffsetLPS "<<IsocenterOffsetLPS<<std::endl;
// std::cout<<"NominalIsocenter IEC"<< m_DRTGeometryMetaInfo->GetProjectionCenter() <<std::endl;
// std::cout<<"CALIBRATION NominalIsocenter_wZcorrectionLPS "<<NominalIsocenter_wZcorrectionLPS<<std::endl;
// std::cout<<"CALIBRATION NominalIsocenterZero_wZcorrectionLPS "<<NominalIsocenterZero_wZcorrectionLPS<<std::endl;
// std::cout<<"CALIBRATION IsocenterOffsetLPS "<<IsocenterOffsetLPS<<std::endl;
ImageType3D::PointType CalibratedIsocenterLPS;
@ -1609,17 +1629,17 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
CalibratedIsocenterLPS[2] = NominalIsocenter_wZcorrectionLPS[2] +
IsocenterOffsetLPS[2];
// std::cout<<"CALIBRATION CalibratedIsocenterLPS "<<CalibratedIsocenterLPS<<std::endl;
// std::cout<<"CALIBRATION CalibratedIsocenterLPS "<<CalibratedIsocenterLPS<<std::endl;
ImageType3D::PointType CalibratedIsocenterZeroLPS;
CalibratedIsocenterZeroLPS[0] = NominalIsocenterZero_wZcorrectionLPS[0] +
IsocenterOffsetLPS[0];
IsocenterOffsetLPS[0] /*+ m_CTMetaInfo->GetSpacing()[0]/2*/;
CalibratedIsocenterZeroLPS[1] = NominalIsocenterZero_wZcorrectionLPS[1] +
IsocenterOffsetLPS[1];
IsocenterOffsetLPS[1] /*+ m_CTMetaInfo->GetSpacing()[1]/2*/;
CalibratedIsocenterZeroLPS[2] = NominalIsocenterZero_wZcorrectionLPS[2] +
IsocenterOffsetLPS[2];
// std::cout<<"CALIBRATION CalibratedIsocenterZeroLPS "<<CalibratedIsocenterZeroLPS<<std::endl;
IsocenterOffsetLPS[2] /*+ m_CTMetaInfo->GetSpacing()[2]/2*/;
// std::cout<<"CALIBRATION CalibratedIsocenterZeroLPS "<<CalibratedIsocenterZeroLPS<<std::endl;
m_DRTImage1MetaInfo->SetProjectionAngleLPS(
this->CalcProjectionAngleLPS(
@ -1903,6 +1923,8 @@ void itkImageProcessor::GetProjectionImages(){
CurrTransform->GetAngleZ()
);
transform2 ->SetCenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero());

8
reg23Topograms/itkDTRrecon/itkImageProcessor.h
View File

@ -152,6 +152,10 @@ public:
/** Panel offsets - panel offsets IEC */
void SetPanelOffsets(double,double);
/** Panel offsets - takes into account patient orientation */
double GetPanelOffsetPGeo(tProjOrientationType ImgPrj);
/** Get projection angles - Gantry angle LPS
* Only meaningful after a 3D Volume is loaded */
double GetPanelOffset1();
@ -311,7 +315,6 @@ private:
toVTKLocalizer2;
InternalImageType::Pointer
image3DIn,
imageDRT1In,
imageDRT2In;
@ -322,9 +325,6 @@ private:
m_PASourceDupli,
m_VolumeSourceDupli;
ChangeInformationFilterType::Pointer
m_3DInputChangeInformationToZero;
/** Apply transform to CT volume.
* This moves the volume based on RT Plan info
* if those are available. */

3
reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h
View File

@ -131,6 +131,9 @@ public:
using ContinuousIndexType = typename Superclass::ContinuousIndexType;
/** Internal floating point comparison accuracy **/
static const double EPSILON;
/** \brief
* Interpolate the image at a point position.
*

797
reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
View File

@ -56,40 +56,45 @@ CIT 6, 89-94 (1998).
namespace itk
{
template <typename TInputImage, typename TCoordRep>
const double gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::EPSILON =
0.00001;
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_PanelOffset = 0.;
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_PanelOffset = 0.;
m_SourcePoint[0] = 0.;
m_SourcePoint[1] = 0.;
m_SourcePoint[2] = 0.;
m_SourcePoint[0] = 0.;
m_SourcePoint[1] = 0.;
m_SourcePoint[2] = 0.;
m_InverseTransform = TransformType::New();
m_InverseTransform->SetComputeZYX(true);
m_InverseTransform = TransformType::New();
m_InverseTransform->SetComputeZYX(true);
m_ComposedTransform = TransformType::New();
m_ComposedTransform->SetComputeZYX(true);
m_ComposedTransform = TransformType::New();
m_ComposedTransform->SetComputeZYX(true);
m_GantryRotTransform = TransformType::New();
m_GantryRotTransform->SetComputeZYX(true);
m_GantryRotTransform->SetIdentity();
m_GantryRotTransform = TransformType::New();
m_GantryRotTransform->SetComputeZYX(true);
m_GantryRotTransform->SetIdentity();
m_CamShiftTransform = TransformType::New();
m_CamShiftTransform->SetComputeZYX(true);
m_CamShiftTransform->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_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;
m_Threshold = 0;
}
@ -97,10 +102,10 @@ template <typename TInputImage, typename TCoordRep>
void
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::PrintSelf(std::ostream & os, Indent indent) const
{
this->Superclass::PrintSelf(os, indent);
this->Superclass::PrintSelf(os, indent);
os << indent << "Threshold: " << m_Threshold << std::endl;
os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
os << indent << "Threshold: " << m_Threshold << std::endl;
os << indent << "Transform: " << m_Transform.GetPointer() << std::endl;
}
@ -108,364 +113,388 @@ template <typename TInputImage, typename TCoordRep>
typename gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::OutputType
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(const PointType & point) const
{
float rayVector[3];
IndexType cIndex;
float rayVector[3];
IndexType cIndex;
PointType drrPixelWorld; // Coordinate of a DRR pixel in the world coordinate system
OutputType pixval;
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;
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();
// 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 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);
}
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 PointReq = point;
PointReq[0] += m_PanelOffset;
PointType PointReq = point;
PointReq[0] += m_PanelOffset;
drrPixelWorld = m_InverseTransform->TransformPoint(PointReq);
drrPixelWorld = m_InverseTransform->TransformPoint(PointReq);
PointType SlidingSourcePoint = m_SourcePoint;
SlidingSourcePoint[0] = 0.;
SlidingSourcePoint[1] = point[1];
SlidingSourcePoint[2] = 0.;
// Get ths input pointers
InputImageConstPointer inputPtr = this->GetInputImage();
PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);
//std::cout<<"SourceWorld: "<<SourceWorld<<std::endl;
typename InputImageType::SizeType sizeCT;
typename InputImageType::RegionType regionCT;
typename InputImageType::SpacingType ctPixelSpacing;
typename InputImageType::PointType ctOrigin;
// 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();
// If Pixel position (in mm) is outside bounds of CT (zero-based)
// assign 0 at the pixel and move on. Ensures regular spacing of resulting
// DRR
float xSizeCT = sizeCT[0] * ctPixelSpacing[0];
float ySizeCT = sizeCT[1] * ctPixelSpacing[1];
float zSizeCT = sizeCT[2] * ctPixelSpacing[2];
float xDrrPix = drrPixelWorld[0];float yDrrPix = drrPixelWorld[1];float zDrrPix = drrPixelWorld[2];
// if(zDrrPix < 0 /*|| yDrrPix < 0 || xDrrPix < 0*/)
// std::cout << drrPixelWorld[0]<<" " <<drrPixelWorld[1]<<" "<<drrPixelWorld[2]<<" / ";
if(xDrrPix>xSizeCT || yDrrPix>ySizeCT || zDrrPix>zSizeCT){
pixval = static_cast<OutputType>(0.0);
return pixval;
}
// calculate the detector position for this pixel center by moving
// 2*m_FocalPointToIsocenterDistance from the source in the pixel
// directions
double p1[3], p2[3];
p1[0] = SourceWorld[0];
p1[1] = SourceWorld[1];
p1[2] = SourceWorld[2];
p2[0] = drrPixelWorld[0];
p2[1] = drrPixelWorld[1];
p2[2] = drrPixelWorld[2];
double P12D = sqrt(
(p2[0]-p1[0]) * (p2[0]-p1[0]) +
(p2[1]-p1[1]) * (p2[1]-p1[1]) +
(p2[2]-p1[2]) * (p2[2]-p1[2])
);
double p12V[3];
p12V[0] = p2[0] - p1[0];
p12V[1] = p2[1] - p1[1];
p12V[2] = p2[2] - p1[2];
double p12v [3];
p12v[0] = p12V[0]/P12D;
p12v[1] = p12V[1]/P12D;
p12v[2] = p12V[2]/P12D;
double pDest [3];
pDest[0] = p1[0] + m_FocalPointToIsocenterDistance * 2 * p12v[0];
pDest[1] = p1[1] + m_FocalPointToIsocenterDistance * 2 * p12v[1];
pDest[2] = p1[2] + m_FocalPointToIsocenterDistance * 2 * p12v[2];
ctPixelSpacing = inputPtr->GetSpacing();
ctOrigin = inputPtr->GetOrigin();
regionCT = inputPtr->GetLargestPossibleRegion();
sizeCT = regionCT.GetSize();
// The following is the Siddon-Jacob fast ray-tracing algorithm
//rayVector[0] = drrPixelWorld[0] - SourceWorld[0];
//rayVector[1] = drrPixelWorld[1] - SourceWorld[1];
PointType SlidingSourcePoint = m_SourcePoint;
SlidingSourcePoint[0] = 0.;
// Simple sliding source model, the source follows the y of the pixel
// requested, which sets the Z position to look in the CT volume.
SlidingSourcePoint[1] = point[1];
SlidingSourcePoint[2] = 0.;
// correct ray to reach detector
rayVector[0] = pDest[0] - SourceWorld[0];
rayVector[1] = pDest[1] - SourceWorld[1];
rayVector[2] = pDest[2] - SourceWorld[2];
PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);
PointType O(3);
O[0] = -ctPixelSpacing[0]/2.;
O[1] = -ctPixelSpacing[1]/2.;
O[2] = -ctPixelSpacing[2]/2.;
// be sure that drr pixel position is inside the volume
float xSizeCT = O[0] + (sizeCT[0]) * ctPixelSpacing[0];
float ySizeCT = O[1] + (sizeCT[1]) * ctPixelSpacing[1];
float zSizeCT = O[2] + (sizeCT[2]) * ctPixelSpacing[2];
float xDrrPix = drrPixelWorld[0];
float yDrrPix = drrPixelWorld[1];
float zDrrPix = drrPixelWorld[2];
if(
xDrrPix<= O[0] || xDrrPix>=xSizeCT ||
yDrrPix<= O[1] || yDrrPix>=ySizeCT ||
zDrrPix<= O[2] || zDrrPix>=zSizeCT
){
pixval = static_cast<OutputType>(0.0);
return pixval;
}
// calculate the detector position for this pixel center by moving
// 2*m_FocalPointToIsocenterDistance from the source in the pixel
// directions
double p1[3], p2[3];
p1[0] = SourceWorld[0];
p1[1] = SourceWorld[1];
p1[2] = SourceWorld[2];
p2[0] = drrPixelWorld[0];
p2[1] = drrPixelWorld[1];
p2[2] = drrPixelWorld[2];
double P12D = sqrt(
(p2[0]-p1[0]) * (p2[0]-p1[0]) +
(p2[1]-p1[1]) * (p2[1]-p1[1]) +
(p2[2]-p1[2]) * (p2[2]-p1[2])
);
double p12V[3];
p12V[0] = p2[0] - p1[0];
p12V[1] = p2[1] - p1[1];
p12V[2] = p2[2] - p1[2];
double p12v [3];
p12v[0] = p12V[0]/P12D;
p12v[1] = p12V[1]/P12D;
p12v[2] = p12V[2]/P12D;
double pDest [3];
pDest[0] = p1[0] + m_FocalPointToIsocenterDistance * 2 * p12v[0];
pDest[1] = p1[1] + m_FocalPointToIsocenterDistance * 2 * p12v[1];
pDest[2] = p1[2] + m_FocalPointToIsocenterDistance * 2 * p12v[2];
/* Calculate the parametric values of the first and the last
// 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 for topogram geometry
// Differs from conventional FPD geometry (above) that has
// the panel at the end of the ray.
rayVector[0] = pDest[0] - SourceWorld[0];
rayVector[1] = pDest[1] - SourceWorld[1];
rayVector[2] = pDest[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))
if (fabs(rayVector[0]) > EPSILON/* != 0*/)
{
/* 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;
// alphaX1 = (0.0 - SourceWorld[0]) / rayVector[0];
// alphaXN = (sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaX1 = ( O[0] - SourceWorld[0]) / rayVector[0];
alphaXN = ( O[0] + sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaXmin = std::min(alphaX1, alphaXN);
alphaXmax = std::max(alphaX1, alphaXN);
}
else
{
/* Current ray front intercepts with z-plane. Update alphaZ. */
alphaCmin = alphaZ;
cIndex[2] = cIndex[2] + kU;
alphaZ = alphaZ + alphaUz;
alphaXmin = -2;
alphaXmax = 2;
}
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 (fabs(rayVector[1]) > EPSILON/* != 0*/)
{
/* 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);
}
alphaY1 = ( O[1] - SourceWorld[1]) / rayVector[1];
alphaYN = ( O[1] + sizeCT[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
// 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 (d12 < minOutputValue)
{
pixval = minOutputValue;
}
else if (d12 > maxOutputValue)
{
pixval = maxOutputValue;
}
else
{
pixval = static_cast<OutputType>(d12);
}
return (pixval);
if (fabs(rayVector[2]) > EPSILON/* != 0*/)
{
// alphaZ1 = (0.0 - SourceWorld[2]) / rayVector[2];
// alphaZN = (sizeCT[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaZ1 = ( O[2] - SourceWorld[2]) / rayVector[2];
alphaZN = ( O[2] + 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] - O[0] ; //ctPixelSpacing[0]/2.;
firstIntersection[1] = SourceWorld[1] + alphaMin * rayVector[1] - O[1] ; //ctPixelSpacing[1]/2.;
firstIntersection[2] = SourceWorld[2] + alphaMin * rayVector[2] - O[2] ; //ctPixelSpacing[2]/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];
// now "correct" the indices (there are special cases):
if (alphaMin == alphaYmin) // j is special
{
if (alphaYmin == alphaY1)
firstIntersectionIndex [1] = 0;
else
// alpha_y_Ny
firstIntersectionIndex [1] = sizeCT[1] - 1;
}
else if (alphaMin == alphaXmin) // i is special
{
if (alphaXmin == alphaX1)
firstIntersectionIndex[0] = 0;
else
// alpha_x_Nx
firstIntersectionIndex[0] = sizeCT[0] - 1;
}
else if (alphaMin == alphaZmin) // k is special
{
if (alphaZmin == alphaZ1)
firstIntersectionIndex [2] = 0;
else
// alpha_z_Nz
firstIntersectionIndex [2] = sizeCT[2] - 1;
}
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 (fabs(rayVector[0]) <= EPSILON/* == 0*/)
{
alphaX = 2;
}
else
{
alphaIntersectionUp[0] =
(O[0] + firstIntersectionIndexUp[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaIntersectionDown[0] =
(O[0] + firstIntersectionIndexDown[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaX = std::max(alphaIntersectionUp[0], alphaIntersectionDown[0]);
}
if (fabs(rayVector[1] ) <= EPSILON/* == 0*/)
{
alphaY = 2;
}
else
{
alphaIntersectionUp[1] =
(O[1] + firstIntersectionIndexUp[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
alphaIntersectionDown[1] =
(O[1] + firstIntersectionIndexDown[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
alphaY = std::max(alphaIntersectionUp[1], alphaIntersectionDown[1]);
}
if (fabs(rayVector[2]) <= EPSILON/* == 0*/)
{
alphaZ = 2;
}
else
{
alphaIntersectionUp[2] =
(O[2] + firstIntersectionIndexUp[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaIntersectionDown[2] =
(O[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 (fabs(rayVector[0]) > EPSILON /*!= 0*/)
{
alphaUx = ctPixelSpacing[0] / fabs(rayVector[0]);
}
else
{
alphaUx = 999;
}
if (fabs(rayVector[1]) > EPSILON /*!= 0*/)
{
alphaUy = ctPixelSpacing[1] / fabs(rayVector[1]);
}
else
{
alphaUy = 999;
}
if (fabs(rayVector[2]) > EPSILON /*!= 0*/)
{
alphaUz = ctPixelSpacing[2] / fabs(rayVector[2]);
}
else
{
alphaUz = 999;
}
/* Calculate voxel index incremental values along the ray path. */
iU = (SourceWorld[0] < drrPixelWorld[0]) ? 1 : -1;
jU = (SourceWorld[1] < drrPixelWorld[1]) ? 1 : -1;
kU = (SourceWorld[2] < drrPixelWorld[2]) ? 1 : -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
const ContinuousIndexType & index) const
{
OutputPointType point;
this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
OutputPointType point;
this->m_Image->TransformContinuousIndexToPhysicalPoint(index, point);
return this->Evaluate(point);
return this->Evaluate(point);
}
@ -473,36 +502,52 @@ template <typename TInputImage, typename TCoordRep>
void
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::ComputeInverseTransform() const
{
m_ComposedTransform->SetIdentity();
m_ComposedTransform->Compose(m_Transform, 0);
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);
typename TransformType::InputPointType isocenter;
isocenter = m_Transform->GetCenter();
// 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);
// //This is problably not necessary since we correct the ray.... to be checked..
// InputImageConstPointer inputPtr = this->GetInputImage();
// 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).)
// std::cout<<inputPtr<<std::endl;
// if(inputPtr != NULL)
// {
// isocenter[0] -= inputPtr->GetSpacing()[0]/2.;
// isocenter[1] -= inputPtr->GetSpacing()[1]/2.;
// isocenter[2] -= inputPtr->GetSpacing()[2]/2.;
// std::cout<<inputPtr->GetSpacing()<<std::endl;
// } else {
// isocenter[0] += 0.5;
// isocenter[1] += 0.5;
// isocenter[2] += 1;
// }
// 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);
m_ComposedTransform->Compose(m_CamRotTransform, 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);
// The overall inverse transform is computed. The inverse transform will be used by the interpolation
// procedure.
m_ComposedTransform->GetInverse(m_InverseTransform);
this->Modified();
// 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();
}
@ -510,8 +555,8 @@ template <typename TInputImage, typename TCoordRep>
void
gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Initialize()
{
this->ComputeInverseTransform();
m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
this->ComputeInverseTransform();
m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
}
} // namespace itk

8
reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.cxx
View File

@ -43,6 +43,7 @@ vtkContourTopogramProjectionFilter::vtkContourTopogramProjectionFilter()
this->dImportOffsetLPS[2] = 0.;
this->dSCD = 0.;
this->dPOffset = 0.;
this->m_RefTransform = nullptr;
this->m_Transform = nullptr;
@ -152,7 +153,7 @@ int vtkContourTopogramProjectionFilter::RequestData(
/** correct x coordinate by the ratio of distance from source
* Basic intercept theorem */
if(pp1[2] != 0.){
pp1[0] = pp1[0] * dSCD / abs(pp1[2]);
pp1[0] = (pp1[0] * dSCD / abs(pp1[2])) + dPOffset;
} else {
vtkDebugMacro("point at SCD distance from panel. Skipping lateral correction for this point.");
}
@ -213,6 +214,11 @@ void vtkContourTopogramProjectionFilter::SetSCD(const double dVal){
dSCD = dVal;
}
void vtkContourTopogramProjectionFilter::SetPanelOffset(const double dVal){
dPOffset = dVal;
}
void vtkContourTopogramProjectionFilter::SetLPStoProjectionLPSGeometryOffset(const double * dP)
{
memcpy(dProjectionLPSOffset,dP,3*sizeof(double));

4
reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.h
View File

@ -40,6 +40,7 @@ public:
void SetLPStoProjectionLPSGeometryOffset (const double * );
void SetSCD(const double dVal);
void SetImportOffsetLPS(const double *);
void SetPanelOffset(const double dVal);
void cleanUpFilter();
@ -60,7 +61,8 @@ private:
double
dProjectionLPSOffset[3],
dImportOffsetLPS[3],
dSCD;
dSCD,
dPOffset;
vtkTransform *m_RefTransform;
vtkTransform *m_Transform;