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

Siddon geometry

Browse Source
This commit is contained in:
Proton local user
2023-06-06 00:04:05 +02:00
parent 1aca4117c6
commit 6f3a5d70b6
2 changed files with 146 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files

82
reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
View File

@ -551,6 +551,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
@ -613,25 +614,49 @@ int itkImageProcessor::fill3DVolumeMeta(
// geometry is fully defined. The origin of the CT image becomes irrelavent.
ImageType3D::PointType pZeroOrigin;
pZeroOrigin[0] = 0.;
pZeroOrigin[1] = 0.;
pZeroOrigin[2] = 0.;
pZeroOrigin[0] = 0;//2312312;//-m_VolumeSourceDupli->GetOutput()->GetOrigin()[0]; //0.;
pZeroOrigin[1] = 0;//23;//-m_VolumeSourceDupli->GetOutput()->GetOrigin()[1]; //0.;
pZeroOrigin[2] = 0;//4444;//-m_VolumeSourceDupli->GetOutput()->GetOrigin()[2]; //0.;
if (m_3DInputChangeInformationToZero) {
m_3DInputChangeInformationToZero = NULL;
m_3DInputChangeInformationToZero = ChangeInformationFilterType::New();
}
m_VolumeSourceDupli->GetOutput()->GetOrigin();
ImageType3D::RegionType lagerReg =
m_VolumeSourceDupli->GetOutput()->GetLargestPossibleRegion();
std::cout<<"CCASASDASASD "<<lagerReg<<std::endl;
ImageType3D::PointType pPoint;
m_VolumeSourceDupli->GetOutput()
->TransformIndexToPhysicalPoint(lagerReg.GetIndex(),pPoint);
std::cout<<"CCASASDASASD "<<pPoint<<std::endl;
//CHECK THIS ZERO.
m_3DInputChangeInformationToZero->SetInput(
m_VolumeSourceDupli->GetOutput());
m_3DInputChangeInformationToZero->SetOutputOrigin(
pZeroOrigin);
m_3DInputChangeInformationToZero->ChangeOriginOn();
m_3DInputChangeInformationToZero->UpdateOutputInformation();
m_3DInputChangeInformationToZero->Update();
image3DIn =
m_3DInputChangeInformationToZero->GetOutput();
lagerReg =
m_3DInputChangeInformationToZero->GetOutput()->GetLargestPossibleRegion();
std::cout<<"CCASASDASASD "<<lagerReg<<std::endl;
m_3DInputChangeInformationToZero->GetOutput()
->TransformIndexToPhysicalPoint(lagerReg.GetIndex(),pPoint);
std::cout<<"CCASASDASASD "<<pPoint<<std::endl;
//m_3DInputChangeInformationToZero->get
return EXIT_SUCCESS;
@ -923,7 +948,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 +967,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 +991,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,6 +1357,8 @@ void itkImageProcessor::InitializeProjector()
transform1->SetCenter(
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() );
std::cout<<"---> CENTER PA " <<
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()<<std::endl;
// 2D Image 1
interpolator1->SetProjectionAngle(
@ -1362,6 +1397,7 @@ void itkImageProcessor::InitializeProjector()
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() );
//transform2->Print(std::cout);
// 2D Image 2
interpolator2->SetProjectionAngle(
dtr * m_DRTImage2MetaInfo->GetProjectionAngleLPS() );
@ -1379,6 +1415,21 @@ void itkImageProcessor::InitializeProjector()
resampleFilter1 = ResampleFilterType::New();
resampleFilter1->SetInput(
image3DIn);
//m_VolumeSourceDupli->GetOutput());
//image3DIn);
ImageType3D::RegionType lagerReg =
m_VolumeSourceDupli->GetOutput()->GetLargestPossibleRegion();
std::cout<<"CCASASDASASD "<<lagerReg<<std::endl;
ImageType3D::PointType pPoint;
m_VolumeSourceDupli->GetOutput()
->TransformIndexToPhysicalPoint(lagerReg.GetIndex(),pPoint);
std::cout<<"CCASASDASASD "<<pPoint<<std::endl;
resampleFilter1->SetDefaultPixelValue(0);
resampleFilter1->SetNumberOfWorkUnits(iNWUnits);
@ -1397,7 +1448,8 @@ void itkImageProcessor::InitializeProjector()
resampleFilter2 = ResampleFilterType::New();
resampleFilter2->SetInput(
image3DIn);
m_VolumeSourceDupli->GetOutput());
//image3DIn);
resampleFilter2->SetDefaultPixelValue(0);
resampleFilter2->SetNumberOfWorkUnits(iNWUnits);
@ -1595,10 +1647,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 +1661,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 +1955,8 @@ void itkImageProcessor::GetProjectionImages(){
CurrTransform->GetAngleZ()
);
transform2 ->SetCenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero());

103
reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
View File

@ -154,13 +154,7 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
drrPixelWorld = m_InverseTransform->TransformPoint(PointReq);
PointType SlidingSourcePoint = m_SourcePoint;
SlidingSourcePoint[0] = 0.;
SlidingSourcePoint[1] = point[1];
SlidingSourcePoint[2] = 0.;
PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);
//std::cout<<"SourceWorld: "<<SourceWorld<<std::endl;
// Get ths input pointers
InputImageConstPointer inputPtr = this->GetInputImage();
@ -175,14 +169,25 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
regionCT = inputPtr->GetLargestPossibleRegion();
sizeCT = regionCT.GetSize();
PointType SlidingSourcePoint = m_SourcePoint;
SlidingSourcePoint[0] = 0.;
SlidingSourcePoint[1] = point[1] /*- ctPixelSpacing[2]/2.*/;
SlidingSourcePoint[2] = 0.;
PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);
//std::cout<<"SourceWorld: "<<SourceWorld<<std::endl;
// 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];
float xSizeCT = -ctPixelSpacing[0]/2. + (sizeCT[0]) * ctPixelSpacing[0];
float ySizeCT = -ctPixelSpacing[1]/2. + (sizeCT[1]) * ctPixelSpacing[1];
float zSizeCT = -ctPixelSpacing[2]/2. + (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]<<" / ";
@ -191,6 +196,15 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
pixval = static_cast<OutputType>(0.0);
return pixval;
}
// std::cout<<point << " " <<PointReq <<" "<<drrPixelWorld<<std::endl;
// std::cout<< ctOrigin <<std::endl;
// std::cout<< m_Transform->GetCenter()<<std::endl;
// std::cout<<sizeCT << " " <<ctPixelSpacing <<" "<<drrPixelWorld<<std::endl;
// exit(1);
// calculate the detector position for this pixel center by moving
// 2*m_FocalPointToIsocenterDistance from the source in the pixel
// directions
@ -241,8 +255,10 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
define the CT volume. */
if (rayVector[0] != 0)
{
alphaX1 = (0.0 - SourceWorld[0]) / rayVector[0];
alphaXN = (sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
// alphaX1 = (0.0 - SourceWorld[0]) / rayVector[0];
// alphaXN = (sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaX1 = ( (-ctPixelSpacing[0]/2.) - SourceWorld[0]) / rayVector[0];
alphaXN = ( (-ctPixelSpacing[0]/2.) + sizeCT[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaXmin = std::min(alphaX1, alphaXN);
alphaXmax = std::max(alphaX1, alphaXN);
}
@ -254,8 +270,10 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
if (rayVector[1] != 0)
{
alphaY1 = (0.0 - SourceWorld[1]) / rayVector[1];
alphaYN = (sizeCT[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
alphaY1 = ( (-ctPixelSpacing[1]/2.) - SourceWorld[1]) / rayVector[1];
alphaYN = ( -ctPixelSpacing[1]/2. + 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);
}
@ -267,8 +285,10 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
if (rayVector[2] != 0)
{
alphaZ1 = (0.0 - SourceWorld[2]) / rayVector[2];
alphaZN = (sizeCT[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
// alphaZ1 = (0.0 - SourceWorld[2]) / rayVector[2];
// alphaZN = (sizeCT[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaZ1 = ( (-ctPixelSpacing[2]/2.) - SourceWorld[2]) / rayVector[2];
alphaZN = ( -ctPixelSpacing[2]/2. + sizeCT[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaZmin = std::min(alphaZ1, alphaZN);
alphaZmax = std::max(alphaZ1, alphaZN);
}
@ -287,15 +307,41 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
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];
firstIntersection[0] = SourceWorld[0] + alphaMin * rayVector[0] - ctPixelSpacing[0]/2.;
firstIntersection[1] = SourceWorld[1] + alphaMin * rayVector[1] - ctPixelSpacing[1]/2.;
firstIntersection[2] = SourceWorld[2] + alphaMin * rayVector[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]);
@ -311,8 +357,8 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
}
else
{
alphaIntersectionUp[0] = (firstIntersectionIndexUp[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaIntersectionDown[0] = (firstIntersectionIndexDown[0] * ctPixelSpacing[0] - SourceWorld[0]) / rayVector[0];
alphaIntersectionUp[0] = (firstIntersectionIndexUp[0] * ctPixelSpacing[0] -(ctPixelSpacing[0]/2.) - SourceWorld[0]) / rayVector[0];
alphaIntersectionDown[0] = (firstIntersectionIndexDown[0] * ctPixelSpacing[0] -(ctPixelSpacing[0]/2.) - SourceWorld[0]) / rayVector[0];
alphaX = std::max(alphaIntersectionUp[0], alphaIntersectionDown[0]);
}
@ -322,8 +368,8 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
}
else
{
alphaIntersectionUp[1] = (firstIntersectionIndexUp[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
alphaIntersectionDown[1] = (firstIntersectionIndexDown[1] * ctPixelSpacing[1] - SourceWorld[1]) / rayVector[1];
alphaIntersectionUp[1] = (firstIntersectionIndexUp[1] * ctPixelSpacing[1] -(ctPixelSpacing[1]/2.) - SourceWorld[1]) / rayVector[1];
alphaIntersectionDown[1] = (firstIntersectionIndexDown[1] * ctPixelSpacing[1] -(ctPixelSpacing[1]/2.) - SourceWorld[1]) / rayVector[1];
alphaY = std::max(alphaIntersectionUp[1], alphaIntersectionDown[1]);
}
@ -333,8 +379,8 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
}
else
{
alphaIntersectionUp[2] = (firstIntersectionIndexUp[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaIntersectionDown[2] = (firstIntersectionIndexDown[2] * ctPixelSpacing[2] - SourceWorld[2]) / rayVector[2];
alphaIntersectionUp[2] = (firstIntersectionIndexUp[2] * ctPixelSpacing[2] -(ctPixelSpacing[2]/2.) - SourceWorld[2]) / rayVector[2];
alphaIntersectionDown[2] = (firstIntersectionIndexDown[2] * ctPixelSpacing[2] -(ctPixelSpacing[2]/2.) - SourceWorld[2]) / rayVector[2];
alphaZ = std::max(alphaIntersectionUp[2], alphaIntersectionDown[2]);
}
@ -478,6 +524,13 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::ComputeInv
typename TransformType::InputPointType isocenter;
isocenter = m_Transform->GetCenter();
// if(this->GetInputImage() != NULL)
//{
// isocenter[0] -= this->GetInputImage()->GetSpacing()[0];
// isocenter[1] -= this->GetInputImage()->GetSpacing()[1];
// isocenter[2] -= this->GetInputImage()->GetSpacing()[2];
//}
// 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.