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

* fixed DRT computation with rotations (6DOF).

Browse Source 
* first step towards input volume alignment form rt plan info
This commit is contained in:
fattori_g
2022-02-11 00:16:51 +01:00
parent 4b8f35b6c4
commit 70498e03b6
3 changed files with 203 additions and 116 deletions
Download Patch File Download Diff File Expand all files Collapse all files

255
reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
View File

@ -46,6 +46,9 @@ gfattori 08.11.2021
#define RT_LNG 0
#define RT_VRT 0
#define TMP_CT3X 4.1305
#define TMP_CT3Y -0.9476
#define TMP_CT3Z 252.7127
namespace itk
{
@ -125,7 +128,7 @@ itkImageProcessor::itkImageProcessor()
image2center[0]=image2center[1]= 0.;
image1Size[0] = image1Size[1] = 512;
image2Size[0] = image2Size[1] = 512;
// image3DCOV[0] = image3DCOV[1] = image3DCOV[2] = 0.;
// image3DCOV[0] = image3DCOV[1] = image3DCOV[2] = 0.;
image1IntensityWindow[0] = image1IntensityWindow[1] = 0.;
image2IntensityWindow[0] = image2IntensityWindow[1] = 0.;
@ -163,6 +166,9 @@ itkImageProcessor::itkImageProcessor()
}
m_3DInputChangeInformation = ChangeInformationFilterInputType::New();
d_CTImportOffset_txiec = CT_IMPORTOFFSET_TXIEC;
d_CTImportOffset_tyiec = CT_IMPORTOFFSET_TYIEC;
d_CTImportOffset_tziec = CT_IMPORTOFFSET_TZIEC;
@ -296,8 +302,8 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
imageSeriesReader3D->Update();
// eImageOrientationType
// currImgOrient;
// eImageOrientationType
// currImgOrient;
/* check patient orientation */
if(fileNames.size()>0){
@ -317,11 +323,11 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::HFS])){
std::cout<<"Image Orientation: HFS"<<std::endl;
// currImgOrient = eImageOrientationType::HFS;
// currImgOrient = eImageOrientationType::HFS;
m_3DPatOrientation = eImageOrientationType::HFS;
} else if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::FFS])){
std::cout<<"Image Orientation: FFS"<<std::endl;
// currImgOrient = eImageOrientationType::FFS;
// currImgOrient = eImageOrientationType::FFS;
m_3DPatOrientation = eImageOrientationType::FFS;
} else {
@ -356,15 +362,66 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
projAngle2 = CalcProjectionAngleLPS(m_3DPatOrientation, projAngle2_IEC);
image3DIn = imageSeriesReader3D->GetOutput();
/* TODO: here is the right place to change origin when required
* after this point, reference info are saved */
ImageType3D::PointType rtCouchOffset;
rtCouchOffset[0]=9.60802;
rtCouchOffset[1]=1333.925;
rtCouchOffset[2]=12.2;
ImageType3D::PointType rtIsocenterLPS;
rtIsocenterLPS[0]= 9.01502;
rtIsocenterLPS[1]= -157.708;
rtIsocenterLPS[2]= -1186.84;
ImageType3D::PointType IEC2DCMMapT;
IEC2DCMMapT[0]= 4.10766744;
IEC2DCMMapT[1]= -2316.44798289;
IEC2DCMMapT[2]= 172.82908123;
ImageType3D::PointType ImportOffset=
this->CalcImportVolumeOffset(rtCouchOffset,
LPStoIEC_Directions,
rtIsocenterLPS,
IEC2DCMMapT);
ImageType3D::PointType Origin;
Origin = imageSeriesReader3D->GetOutput()->GetOrigin();
std::cout<< "OldOrigin: " <<Origin<<std::endl;
Origin = Origin - ImportOffset ;
std::cout<< "NewOrigin: " <<Origin<<std::endl;
// exit(1);
// ImageType3D::PointType NewOrigin;
// ImageType3D::PointType Offset;
// Offset[0] = -TMP_CT3X;
// Offset[1] = -TMP_CT3Y;
// Offset[2] = -TMP_CT3Z;
// NewOrigin = imageSeriesReader3D->GetOutput()->GetOrigin();
// std::cout<< "OldOrigin: " <<NewOrigin<<std::endl;
// NewOrigin = NewOrigin - Offset;
m_3DInputChangeInformation->SetInput(imageSeriesReader3D->GetOutput());
m_3DInputChangeInformation->SetOutputOrigin(Origin ) ;//NewOrigin);
m_3DInputChangeInformation->ChangeOriginOn();
m_3DInputChangeInformation->UpdateOutputInformation();
m_3DInputChangeInformation->Update();
std::cout<< "NewOrigin: " <<m_3DInputChangeInformation->GetOutput()->GetOrigin() <<std::endl;
/* End of origin change */
image3DIn = m_3DInputChangeInformation->GetOutput();
/* before changing origin for projection, let's save the center of volume */
using ImageRegionType3D = ImageType3D::RegionType;
using SizeType3D = ImageRegionType3D::SizeType;
@ -372,7 +429,7 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
ImageRegionType3D region3D = image3DIn->GetBufferedRegion();
SizeType3D size3D = region3D.GetSize();
// ImageType3D::PointType origin3D = image3DIn->GetOrigin();
// ImageType3D::PointType origin3D = image3DIn->GetOrigin();
const itk::Vector<double, 3> resolution3D = image3DIn->GetSpacing();
ImageDirectionType3D imagDirection = image3DIn->GetDirection();
@ -398,27 +455,28 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
/* Calculate center of projection in patient coordinates */
/* center of projection in an IEC reference at Patient Origin */
ImageType3D::PointType ProjectionCenterFixedAxes;
ProjectionCenterFixedAxes[0] = 0.;
ProjectionCenterFixedAxes[1] = 0.;
ProjectionCenterFixedAxes[2] = CT_DEFAULT_TABLEHEIGHT;
ImageType3D::PointType ProjectionCenterFixedAxes;
ProjectionCenterFixedAxes[0] = 0.;
ProjectionCenterFixedAxes[1] = 0.;
ProjectionCenterFixedAxes[2] = CT_DEFAULT_TABLEHEIGHT;
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions = LPStoIEC_Directions.GetTranspose();
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions = LPStoIEC_Directions.GetTranspose();
m_3DProjectionOriginLPS =
IECtoLPS_Directions * ProjectionCenterFixedAxes ;
/* in longitudinal direction (Sup-Inf), we put it in the
m_3DProjectionOriginLPS =
IECtoLPS_Directions * ProjectionCenterFixedAxes ;
/* in longitudinal direction (Sup-Inf), we put it in the
* middle of the volume */
m_3DProjectionOriginLPS[2] = image3DCOVLPS[2];
m_3DProjectionOriginLPS[2] = image3DCOVLPS[2];
/* calculate also the same center with zero origin
/* calculate also the same center with zero origin
* This is required to bring back the DRT into LPS */
m_3DProjectionOriginLPSZero =
m_3DProjectionOriginLPS - m_3DOriginLPS;
m_3DProjectionOriginLPSZero =
m_3DProjectionOriginLPS - m_3DOriginLPS;
/* We should know everything we need to calculate the
* origin and direction of projection images */
if(true){
@ -705,22 +763,6 @@ int itkImageProcessor::load2D(const char * pcFName){
InternalImageType::Pointer MyLocalizerImage =
caster2DInput->GetOutput();
// ChangeInformationFilterType::Pointer filter = ChangeInformationFilterType::New();
// filter->SetInput(MyLocalizerImage);
// try
// {
// filter->UpdateOutputInformation();
// }
// catch (itk::ExceptionObject & error)
// {
// std::cerr << "Error: " << error << std::endl;
// return EXIT_FAILURE;
// }
// filter->Update();
double* m_ImageIntensity;
FlipFilterType::Pointer m_FlipFilter;
DuplicatorType::Pointer m_Duplicator;
@ -731,7 +773,7 @@ int itkImageProcessor::load2D(const char * pcFName){
m_ImageIntensity = image1IntensityWindow;
m_Duplicator = m_PASourceDupli;
// m_ImLoaded = &image2D1Loaded;
// m_ImLoaded = &image2D1Loaded;
break;
@ -739,7 +781,7 @@ int itkImageProcessor::load2D(const char * pcFName){
m_ImageIntensity = image2IntensityWindow;
m_Duplicator = m_LATSourceDupli;
// m_ImLoaded = &image2D2Loaded;
// m_ImLoaded = &image2D2Loaded;
break;
@ -788,6 +830,7 @@ void itkImageProcessor::InitializeProjector()
{
caster3D->Update();
transform->SetComputeZYX(true);
TransformType::OutputVectorType translation;
translation[0] = TZero[0];
translation[1] = TZero[1];
@ -802,7 +845,7 @@ void itkImageProcessor::InitializeProjector()
// std::cout<<" -------- PROJECTOR GEOMETRY --------"<<std::endl;
std::cout<<"Projector Isocenter "<<isocenter3D[0]<<" "<<isocenter3D[1]<<" "<<isocenter3D[2] <<std::endl;
// std::cout<<origin3D[0]<<" "<<resolution3D[0]<<" "<<static_cast<double>(size3D[0])<<" "<<static_cast<double>(size3D[0]) / 2.0<<std::endl;
// std::cout<<origin3D[0]<<" "<<resolution3D[0]<<" "<<static_cast<double>(size3D[0])<<" "<<static_cast<double>(size3D[0]) / 2.0<<std::endl;
// std::cout<<"Projector scd "<<scd <<std::endl;
// std::cout<<" -------- -------- --------"<<std::endl;
@ -926,6 +969,8 @@ void itkImageProcessor::GetProjectionImages(){
resampleFilter1->SetInput(caster3D->GetOutput()); // Link the 3D volume.
resampleFilter1->SetDefaultPixelValue(0);
resampleFilter1->SetNumberOfWorkUnits(1);
// The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance
// have been set before registration. Here we only need to replace the initial
// transform with the final transform.
@ -939,36 +984,27 @@ void itkImageProcessor::GetProjectionImages(){
double spacing[Dimension];
double o2Dx, o2Dy, origin[Dimension];
if(false){
// if(image2D1Loaded){
// resampleFilter1->SetSize(m_Input2DImage_PA ->GetLargestPossibleRegion().GetSize());
// resampleFilter1->SetOutputOrigin(m_Input2DImage_PA ->GetOrigin());
// resampleFilter1->SetOutputSpacing(m_Input2DImage_PA ->GetSpacing());
// std::cout<<"Size "<<caster2D1->GetOutput()->GetLargestPossibleRegion().GetSize()<<std::endl;
// std::cout<<"Origin "<<caster2D1->GetOutput()->GetOrigin()<<std::endl;
// std::cout<<"Spacing "<<caster2D1->GetOutput()->GetSpacing()<<std::endl;
} else {
size[0] = image1Size[0]; // number of pixels along X of the 2D DRR image
size[1] = image1Size[1]; // number of pixels along X of the 2D DRR image
size[2] = 1; // only one slice
spacing[0] = image1res[0]; // pixel spacing along X of the 2D DRR image [mm]
spacing[1] = image1res[1]; // pixel spacing along Y of the 2D DRR image [mm]
spacing[2] = 1.0; // slice thickness of the 2D DRR image [mm]
size[0] = image1Size[0]; // number of pixels along X of the 2D DRR image
size[1] = image1Size[1]; // number of pixels along X of the 2D DRR image
size[2] = 1; // only one slice
spacing[0] = image1res[0]; // pixel spacing along X of the 2D DRR image [mm]
spacing[1] = image1res[1]; // pixel spacing along Y of the 2D DRR image [mm]
spacing[2] = 1.0; // slice thickness of the 2D DRR image [mm]
o2Dx = ((double)image1Size[0] - 1.) / 2.;
o2Dy = ((double)image1Size[1] - 1.) / 2.;
// Compute the origin (in mm) of the 2D Image
origin[0] = -image1res[0] * o2Dx;
origin[1] = -image1res[1] * o2Dy;
origin[2] = -scd;
resampleFilter1->SetSize(size);
resampleFilter1->SetOutputSpacing(spacing);
resampleFilter1->SetOutputOrigin(origin);
o2Dx = ((double)image1Size[0] - 1.) / 2.;
o2Dy = ((double)image1Size[1] - 1.) / 2.;
// Compute the origin (in mm) of the 2D Image
origin[0] = -image1res[0] * o2Dx;
origin[1] = -image1res[1] * o2Dy;
origin[2] = -scd;
resampleFilter1->SetSize(size);
resampleFilter1->SetOutputSpacing(spacing);
resampleFilter1->SetOutputOrigin(origin);
}
// Do the same thing for the output image 2.
@ -976,6 +1012,10 @@ void itkImageProcessor::GetProjectionImages(){
resampleFilter2->SetInput(caster3D->GetOutput());
resampleFilter2->SetDefaultPixelValue(0);
resampleFilter1->SetNumberOfWorkUnits(1);
// The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance
// have been set before registration. Here we only need to replace the initial
// transform with the final transform.
@ -983,33 +1023,27 @@ void itkImageProcessor::GetProjectionImages(){
interpolator2->Initialize();
resampleFilter2->SetInterpolator(interpolator2);
if(false){
// if(image2D2Loaded){
// resampleFilter2->SetSize(m_Input2DImage_LAT->GetLargestPossibleRegion().GetSize());
// resampleFilter2->SetOutputOrigin(m_Input2DImage_LAT ->GetOrigin());
// resampleFilter2->SetOutputSpacing(m_Input2DImage_LAT ->GetSpacing());
} else {
size[0] = image2Size[0]; // number of pixels along X of the 2D DRR image
size[1] = image2Size[1]; // number of pixels along X of the 2D DRR image
size[2] = 1; // only one slice
spacing[0] = image2res[0]; // pixel spacing along X of the 2D DRR image [mm]
spacing[1] = image2res[1]; // pixel spacing along Y of the 2D DRR image [mm]
spacing[2] = 1.0; // slice thickness of the 2D DRR image [mm]
size[0] = image2Size[0]; // number of pixels along X of the 2D DRR image
size[1] = image2Size[1]; // number of pixels along X of the 2D DRR image
size[2] = 1; // only one slice
spacing[0] = image2res[0]; // pixel spacing along X of the 2D DRR image [mm]
spacing[1] = image2res[1]; // pixel spacing along Y of the 2D DRR image [mm]
spacing[2] = 1.0; // slice thickness of the 2D DRR image [mm]
o2Dx = ((double)image2Size[0] - 1.) / 2.;
o2Dy = ((double)image2Size[1] - 1.) / 2.;
// Compute the origin (in mm) of the 2D Image
origin[0] = -image2res[0] * o2Dx;
origin[1] = -image2res[1] * o2Dy;
origin[2] = -scd;
o2Dx = ((double)image2Size[0] - 1.) / 2.;
o2Dy = ((double)image2Size[1] - 1.) / 2.;
// Compute the origin (in mm) of the 2D Image
origin[0] = -image2res[0] * o2Dx;
origin[1] = -image2res[1] * o2Dy;
origin[2] = -scd;
resampleFilter2->SetSize(size);
resampleFilter2->SetOutputSpacing(spacing);
resampleFilter2->SetOutputOrigin(origin);
resampleFilter2->SetSize(size);
resampleFilter2->SetOutputSpacing(spacing);
resampleFilter2->SetOutputOrigin(origin);
}
resampleFilter1->Update();
resampleFilter2->Update();
@ -1033,9 +1067,7 @@ void itkImageProcessor::GetProjectionImages(){
/* First of all we need the center of the Projection images in its reference frame */
ImageType3D::PointType image3DOrigin = m_3DProjectionOriginLPS;
// image3DOrigin[0]=image3DCOV[0];
// image3DOrigin[1]=image3DCOV[1];
// image3DOrigin[2]=image3DCOV[2];
resampleFilter1->Update();
@ -1055,8 +1087,8 @@ void itkImageProcessor::GetProjectionImages(){
// InternalImageType::DirectionType LATtoHFS =
// this->CalcDRTImageDirections(projAngle2);
// InternalImageType::DirectionType LATtoHFS =
// this->CalcDRTImageDirections(projAngle2);
/* Again.. */
@ -1087,6 +1119,36 @@ void itkImageProcessor::GetProjectionImages(){
}
itkImageProcessor::ImageType3D::PointType
itkImageProcessor::CalcImportVolumeOffset(
ImageType3D::PointType rtCouchOffset,
InternalImageType::DirectionType VolumeLPStoIEC,
ImageType3D::PointType rtIsocenterLPS,
ImageType3D::PointType IEC2DCMMapT){
InternalImageType::DirectionType VolumeIECtoLPS;
VolumeIECtoLPS = VolumeLPStoIEC.GetTranspose();
ImageType3D::PointType IsoSupport_IECPos;
IsoSupport_IECPos[0] = rtCouchOffset[0] + IEC2DCMMapT[0];
IsoSupport_IECPos[1] = rtCouchOffset[1] + IEC2DCMMapT[1];
IsoSupport_IECPos[2] = rtCouchOffset[2] + IEC2DCMMapT[2];
ImageType3D::PointType IsoSupport_LPSPos =
VolumeIECtoLPS * IsoSupport_IECPos;
ImageType3D::PointType Offset =
rtIsocenterLPS - IsoSupport_LPSPos;
return
Offset;
}
void itkImageProcessor::Write2DImages(){
@ -1162,7 +1224,7 @@ vtkImageData* itkImageProcessor::GetLocalizer1VTK()
toVTKLocalizer1->SetInput(rescaler1->GetOutput());
toVTKLocalizer1->Update();
if(true) {
if(false) {
using ImageRegionType3D = ImageType3D::RegionType;
using SizeType3D = ImageRegionType3D::SizeType;
using ImageDirectionType3D = ImageType3D::DirectionType;
@ -1276,7 +1338,6 @@ vtkImageData* itkImageProcessor::GetProjection1VTK()
rescaler1->SetOutputMinimum(0);
rescaler1->SetOutputMaximum(255);
rescaler1->SetInput( imageDRT1In );
// flipFilter1->GetOutput());
rescaler1->Update();
@ -1333,9 +1394,10 @@ vtkImageData* itkImageProcessor::GetProjection2VTK()
rescaler2->SetOutputMinimum(0);
rescaler2->SetOutputMaximum(255);
rescaler2->SetInput( imageDRT2In );
// flipFilter2->GetOutput());
rescaler2->Update();
toVTK2D2->SetInput(rescaler2->GetOutput());
toVTK2D2->Update();
// rescaler2->Print(std::cout);
@ -1367,8 +1429,7 @@ vtkImageData* itkImageProcessor::GetProjection2VTK()
// std::cout<<"Last Voxel position: "<<LastVoxelPosition<<std::endl;
toVTK2D2->SetInput(rescaler2->GetOutput());
toVTK2D2->Update();
// double* dBounds = toVTK2D2->GetOutput()->GetBounds();

18
reg23Topograms/itkDTRrecon/itkImageProcessor.h
View File

@ -105,7 +105,6 @@ public:
itkTypeMacro(itkImageProcessor, Object);
/** Input data load methods*/
//void load3DVolume(const char *);
int load3DSerie(const char* );
int load2D(const char *);
@ -251,6 +250,10 @@ private:
DuplicatorType::Pointer m_LATSourceDupli;
DuplicatorType::Pointer m_PASourceDupli;
ChangeInformationFilterInputType::Pointer
m_3DInputChangeInformation;
tPatOrientation
m_3DPatOrientation;
@ -278,13 +281,12 @@ private:
tPatOrientation pOrient,
double pAngleIEC);
// ImageType3D::PointType
// CalcDataCollectionCenterLPS_ZeroOrigin(
// tPatOrientation pOrient,
// ImageType3D::PointType m_CollectionCenterIEC,
// ImageType3D::PointType m_CTOriginLPS);
ImageType3D::PointType
CalcImportVolumeOffset(
ImageType3D::PointType rtCouchOffset,
InternalImageType::DirectionType VolumeLPStoIEC,
ImageType3D::PointType rtIsocenterLPS,
ImageType3D::PointType IEC2DCMMapT);
double image1res[2], image2res[2];

46
reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
View File

@ -148,9 +148,18 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
// m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
}
//PointType m_SourcePointSliding = m_SourcePoint;
//m_SourcePointSliding[2] = point[2];
PointType SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
//std::cout<<"drrPixelWorld: "<< point <<std::endl;
drrPixelWorld = m_InverseTransform->TransformPoint(point);
//std::cout<<"drrPixelWorld: "<< drrPixelWorld <<std::endl;
//std::cout<<"m_SourcePoint : "<<m_SourcePoint <<std::endl;
PointType SlidingSourcePoint = m_SourcePoint;
SlidingSourcePoint[1] = point[1];
//PointType SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);
//std::cout<<"SourceWorld: "<<SourceWorld<<std::endl;
// Get ths input pointers
InputImageConstPointer inputPtr = this->GetInputImage();
@ -165,34 +174,44 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
regionCT = inputPtr->GetLargestPossibleRegion();
sizeCT = regionCT.GetSize();
drrPixelWorld = m_InverseTransform->TransformPoint(point);
//slide source along z axis...
SourceWorld[2] = drrPixelWorld[2];
//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];
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]) );
double p12V[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 [2];
double p12v [3];
p12v[0] = p12V[0]/P12D;
p12v[1] = p12V[1]/P12D;
p12v[2] = p12V[2]/P12D;
double pDest [2];
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];
// The following is the Siddon-Jacob fast ray-tracing algorithm
@ -202,8 +221,13 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
// correct ray to reach detector
rayVector[0] = pDest[0] - SourceWorld[0];
rayVector[1] = pDest[1] - SourceWorld[1];
rayVector[2] = drrPixelWorld[2] - SourceWorld[2];
rayVector[2] = pDest[2] - SourceWorld[2];
//rayVector[2] = drrPixelWorld[2] - SourceWorld[2];
//std::cout<< "rayVector "<<rayVector <<std::endl;
//return 0;
/* Calculate the parametric values of the first and the last
intersection points of the ray with the X, Y, and Z-planes that