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

Merge branch 'ScoutRTRelease' into 'master'

Browse Source 
Merging First Release into Master

See merge request cpt_bioeng/drt!37
This commit is contained in:
via_r
2023-08-30 10:35:42 +02:00
parent 14f4aca6d1 15a5154ec5
commit 632303efcc
18 changed files with 1371 additions and 593 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
reg23Topograms/itkDTRrecon/CMakeLists.txt
View File

@ -20,6 +20,7 @@ SET(SRCS
vtkContourTopogramProjectionFilter.cxx
DRTMetaInformation.cpp
itkReg23.cpp
itkReg23MetaInformation.cpp
)
SET(HDR
@ -31,6 +32,7 @@ SET(HDR
vtkContourTopogramProjectionFilter.h
DRTMetaInformation.h
itkReg23.h
itkReg23MetaInformation.h
)
ADD_LIBRARY(${LIB_NAME} ${SRCS} ${HDR})

3
reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp
View File

@ -94,6 +94,8 @@ DRTImageMetaInformation(){
this->m_OriginLPS.Fill(0.);
this->m_PanelOffsetPGeo = 0.;
}
void
@ -549,7 +551,6 @@ R23MetaInformation (){
m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN;
m_OptimizerType = tOptimizerTypeEnum::POWELL;
m_MetricType = tMetricTypeEnum::NCC;
m_MaxIterations = 6;
}

61
reg23Topograms/itkDTRrecon/DRTMetaInformation.h
View File

@ -24,27 +24,34 @@ typedef enum eImageOrientationType{
FFS = 2
} tPatOrientation;
typedef enum eHandlingRotImpTransform {
ALWAYS_USE = 0,
NEVER_USE =1,
STATION_DEPENDENT=2
} tHandlingRotImpTransform;
typedef enum eDegreeOfFreedomType {
UNKNOWN = 0,
THREE_DOF,
SIX_DOF,
X_ONLY,
Y_ONLY,
Z_ONLY,
ROTATION_ONLY,
OTHER
THREE_DOF =1,
SIX_DOF =2 ,
X_ONLY =3 ,
Y_ONLY =4,
Z_ONLY =5 ,
ROTATION_ONLY =6,
OTHER =7
} tDegreeOfFreedomEnum;
typedef enum eOptimizerType{
POWELL,
AMOEBA,
EXHAUSTIVE
POWELL = 0,
AMOEBA = 1,
EXHAUSTIVE =2
} tOptimizerTypeEnum;
typedef enum eMetricType{
NCC,
MI
NCC = 0,
MI = 1
} tMetricTypeEnum;
@ -66,7 +73,8 @@ inline int GetNumberOfDOF(eDegreeOfFreedomType dof)
default:
return 0;
}
}
return 0;
};
class DegreeOfFreedom : public itk::Object {
@ -197,6 +205,9 @@ public:
itkSetMacro(ProjectionOriginLPSZero,PointType);
itkGetMacro(ProjectionOriginLPSZero,PointType);
itkSetMacro(PanelOffsetPGeo,double);
itkGetMacro(PanelOffsetPGeo,double);
itkSetMacro(ProjectionAngleLPS,double);
itkGetMacro(ProjectionAngleLPS,double);
@ -229,7 +240,8 @@ protected:
double
m_ProjectionAngleLPS,
m_SCD;
m_SCD,
m_PanelOffsetPGeo;
DirectionType
m_ImageDirectionsLPS;
@ -424,6 +436,11 @@ public:
itkSetMacro(UseRotationsForHiddenTransform, bool);
itkGetMacro(UseRotationsForHiddenTransform, bool);
itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum);
itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum);
itkSetEnumMacro(HandleRotationImpOffset, tHandlingRotImpTransform);
itkGetEnumMacro(HandleRotationImpOffset, tHandlingRotImpTransform);
protected:
@ -458,7 +475,7 @@ protected:
PointType
/* center of projection in an IEC reference at
* Patient Origin of fixed images. Positioning scanner */
* Patient Origin of fixed images. Positionin scanner */
m_ProjectionCenter,
m_ProjectionCenterOffset1,
m_ProjectionCenterOffset2;
@ -467,6 +484,12 @@ protected:
bool
m_UseRotationsForHiddenTransform;
tDegreeOfFreedomEnum
m_DegreeOfFreedom;
tHandlingRotImpTransform
m_HandleRotationImpOffset;
/** Default Constructor **/
DRTProjectionGeometryImageMetaInformation ();
/** Default Destructor **/
@ -622,8 +645,7 @@ public:
itkSetEnumMacro(MetricType, tMetricTypeEnum);
itkGetEnumMacro(MetricType, tMetricTypeEnum);
itkSetEnumMacro(MaxIterations, int);
itkGetEnumMacro(MaxIterations, int);
protected:
@ -637,9 +659,6 @@ protected:
tMetricTypeEnum
m_MetricType;
int
m_MaxIterations;
/** Default Constructor **/
R23MetaInformation ();
/** Default Destructor **/
@ -712,6 +731,6 @@ private:
void operator=(const Self&);
};
}
};
#endif

2
reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.h
View File

@ -95,6 +95,7 @@ public:
* typically results in narrower valleys in the cost fucntion.
* Default value is false. */
itkSetMacro(SubtractMean, bool);
itkSetMacro(NumberOfHistogramBins, int);
itkGetConstReferenceMacro(SubtractMean, bool);
itkBooleanMacro(SubtractMean);
@ -106,6 +107,7 @@ protected:
private:
bool m_SubtractMean;
int m_NumberOfHistogramBins;
};
} // end namespace itk

6
reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx
View File

@ -36,6 +36,7 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage,
TMovingImage>::MutualInformationTwoImageToOneImageMetric()
{
m_SubtractMean = false;
m_NumberOfHistogramBins = 50;
}
template <typename TFixedImage, typename TMovingImage>
@ -147,7 +148,8 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue()
//auto numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
// since we have a ROI, then we should not set allPixels to TRUE.
//metric1->UseAllPixelsOn();
metric1->SetNumberOfHistogramBins(50);
//std::cout << "m_NumberOfHistogramBins " << m_NumberOfHistogramBins << std::endl;
metric1->SetNumberOfHistogramBins(m_NumberOfHistogramBins);
// InternalImageType::IndexType pIndex;
// pIndex[0]=200;
@ -183,7 +185,7 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue()
//numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
//metric2->SetNumberOfSpatialSamples(numberOfSamples);
//metric2->UseAllPixelsOn();
metric2->SetNumberOfHistogramBins(50);
metric2->SetNumberOfHistogramBins(m_NumberOfHistogramBins);
metric2->SetFixedImage(fixedImage2);
metric2->SetMovingImage(this->m_Filter2->GetOutput());
metric2->Initialize();

400
reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
View File

@ -1,4 +1,4 @@

/*
@ -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;
@ -98,6 +92,18 @@ itkImageProcessor::itkImageProcessor()
m_r23MetaInfo = NULL;
m_r23MetaInfo = R23MetaInformation::New();
m_PowellOptimizerMetaInfo = NULL;
m_PowellOptimizerMetaInfo = PowellOptimizerMetaInformation::New();
m_AmoebaOptimizerMetaInfo = NULL;
m_AmoebaOptimizerMetaInfo = AmoebaOptimizerMetaInformation::New();
m_MIMetricMetaInfo = NULL;
m_MIMetricMetaInfo = MIMetricMetaInformation::New();
m_NCCMetricMetaInfo = NULL;
m_NCCMetricMetaInfo = NCCMetricMetaInformation::New();
/* Initialise the projection geoemtry with defaults */
m_DRTGeometryMetaInfo = NULL;
m_DRTGeometryMetaInfo
@ -108,7 +114,7 @@ itkImageProcessor::itkImageProcessor()
m_DRTGeometryMetaInfo->SetProjectionAngle2IEC(90.);
m_DRTGeometryMetaInfo->SetIntensityThreshold(-300.);
ImageType3D::PointType
Point3D;
Point3D(3);
Point3D[0]=0.;
Point3D[1]=0.;
Point3D[2]=175.;
@ -122,7 +128,7 @@ itkImageProcessor::itkImageProcessor()
m_DRTGeometryMetaInfo->SetDRT1Size(ImageSize);
m_DRTGeometryMetaInfo->SetDRT2Size(ImageSize);
ImageType3D::SpacingType
ImageRes;
ImageRes(3);
ImageRes[0] = 1.;
ImageRes[1] = 1.;
ImageRes[2] = 1.;
@ -136,16 +142,15 @@ itkImageProcessor::itkImageProcessor()
m_DRTGeometryMetaInfo->SetProjectionCenterOffset1(Point3D);
m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(Point3D);
m_DRTGeometryMetaInfo->SetDegreeOfFreedom(eDegreeOfFreedomType::THREE_DOF);
m_DRTGeometryMetaInfo->SetHandleRotationImpOffset(eHandlingRotImpTransform::ALWAYS_USE);
optimizerObserver = CommandIterationUpdate::New();
ExhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New();
// TOOOOODOOOOO optimizerObserver->SetProcess(registration);
m_R23 = itkReg23::New();
m_R23->SetOptimizerObserver(optimizerObserver);
}
itkImageProcessor::~itkImageProcessor()
@ -216,6 +221,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);
@ -388,6 +425,8 @@ int itkImageProcessor::unload3DVolumeAndMeta(){
m_TransformMetaInfo = TransformMetaInformation::New();
this->ResetROIRegions();
return 1;
}
@ -466,12 +505,78 @@ int itkImageProcessor::load3DSerieFromFiles( std::vector<std::string> v_fnames){
return -1;
}
////////////////////////*****************//////////////////////////
CastFilterType3D::Pointer caster3D =
CastFilterType3D::New();
caster3D->SetInput(imageSeriesReader3D->GetOutput());
caster3D->Update();
if(m_ResamplepCT == false){
caster3D->SetInput(imageSeriesReader3D->GetOutput());
caster3D->Update();
}else{
std::cout << "Resampling pCT @ spacing in Z " << m_SamplingLNG << std::endl;
LinearInterpolatorType::Pointer
linearInterpolator = LinearInterpolatorType::New();
// Set identity transform
TransformType::Pointer
transform = TransformType::New();
transform->SetIdentity();
const typename ImageType3D::SpacingType& inputSpacing =
imageSeriesReader3D->GetOutput()->GetSpacing();
const typename ImageType3D::RegionType& inputRegion =
imageSeriesReader3D->GetOutput()->GetLargestPossibleRegion();
const typename ImageType3D::SizeType& inputSize =
inputRegion.GetSize();
// Inputs Info
// std::cout << "inputSpacing " << std::endl;
// std::cout << inputSpacing << std::endl;
// std::cout << "inputRegion " << std::endl;
// std::cout << inputRegion << std::endl;
// std::cout << "inputSize " << std::endl;
// std::cout << inputSize << std::endl;
ImageType3D::SpacingType outputSpacing;
outputSpacing[0] = inputSpacing[0];
outputSpacing[1] = inputSpacing[1];
outputSpacing[2] = m_SamplingLNG;
ImageType3D::SizeType outputSize;
// typedef typename ImageType3D::SizeType::SizeValueType SizeValueType;
outputSize[0] = static_cast<SizeValueType>(
inputSize[0] * inputSpacing[0] / outputSpacing[0] + .5);
outputSize[1] = static_cast<SizeValueType>(
inputSize[1] * inputSpacing[1] / outputSpacing[1] + .5);
outputSize[2] = static_cast<SizeValueType>(
inputSize[2] * inputSpacing[2] / outputSpacing[2] + .5);
ResampleInputFilterType::Pointer
resampler = ResampleInputFilterType::New();
resampler->SetInput( imageSeriesReader3D->GetOutput() );
resampler->SetTransform( transform );
resampler->SetInterpolator( linearInterpolator );
resampler->SetDefaultPixelValue(-1024);
resampler->SetOutputOrigin( imageSeriesReader3D->GetOutput()->GetOrigin() );
resampler->SetOutputSpacing( outputSpacing );
resampler->SetOutputDirection( imageSeriesReader3D->GetOutput()->GetDirection() );
resampler->SetSize( outputSize );
resampler->Update();
caster3D->SetInput(resampler->GetOutput());
caster3D->Update();
}
////////////////////////*****************//////////////////////////
if (m_VolumeSourceDupli == NULL)
std::cout << "NEVER HERE m_VolumeSourceDupli" << std::endl;
@ -535,6 +640,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
@ -596,25 +702,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;
@ -658,6 +747,31 @@ const std::vector <double> itkImageProcessor::GetLPStoProjectionGeoLPSOffset()
vOffset;
}
double
itkImageProcessor::CalcPanelOffsetPGeo(
tPatOrientation pOrient,
double dPOffset){
double dVal = dPOffset;//0.;
switch (pOrient) {
case tPatOrientation :: HFS:
break;
case tPatOrientation :: FFS:
dVal *= -1;
break;
default:
break;
}
return
dVal;
}
double
itkImageProcessor::CalcProjectionAngleLPS(
tPatOrientation pOrient,
@ -728,8 +842,6 @@ void itkImageProcessor::SetProjectionAngleOffsetIEC(double dOff1, double dOff2)
{
m_DRTGeometryMetaInfo->SetProjectionAngle1OffsetIEC(dOff1);
m_DRTGeometryMetaInfo->SetProjectionAngle2OffsetIEC(dOff2);
std::cout << "SetProjectionAngleOffsetIEC " << dOff1 << " " << dOff2 << std::endl;
}
@ -884,7 +996,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;
@ -903,11 +1015,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 =
@ -921,6 +1039,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;
@ -1051,6 +1171,7 @@ double itkImageProcessor::GetLocalizerDisplayWindowLevel(int iImg)
m_TImageMeta->GetWLLevel();
}
return 0.;
}
double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg)
@ -1081,6 +1202,7 @@ double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg)
m_TImageMeta->GetWLWindow();
}
return 0.;
}
@ -1112,18 +1234,20 @@ itkImageProcessor::GetFinalR23Parameters(){
break;
case SIX_DOF:
pPars[0] = currentPosition[0]
- m_TransformMetaInfo->GetHiddenRotations()[0];
pPars[1] = currentPosition[1]
- m_TransformMetaInfo->GetHiddenRotations()[1];
pPars[2] = currentPosition[2]
- m_TransformMetaInfo->GetHiddenRotations()[2];
pPars[3] = currentPosition[3]
pPars[3] = currentPosition[0]
- m_TransformMetaInfo->GetHiddenTranslations()[0];
pPars[4] = currentPosition[4]
pPars[4] = currentPosition[1]
- m_TransformMetaInfo->GetHiddenTranslations()[1];
pPars[5] = currentPosition[5]
pPars[5] = currentPosition[2]
- m_TransformMetaInfo->GetHiddenTranslations()[2];
pPars[0] = currentPosition[3]
- m_TransformMetaInfo->GetHiddenRotations()[0];
pPars[1] = currentPosition[4]
- m_TransformMetaInfo->GetHiddenRotations()[1];
pPars[2] = currentPosition[5]
- m_TransformMetaInfo->GetHiddenRotations()[2];
break;
default:
@ -1137,6 +1261,27 @@ itkImageProcessor::GetFinalR23Parameters(){
}
void itkImageProcessor::SetOptimizer(tOptimizerTypeEnum eOpti){
if(eOpti != tOptimizerTypeEnum::AMOEBA &&
eOpti != tOptimizerTypeEnum::POWELL){
itkExceptionMacro(<< "Unkown optimzer : " << eOpti);
return;
}
std::cout<<"Setting Optimizer: "<<eOpti<<std::endl;
m_r23MetaInfo->SetOptimizerType(eOpti);
}
void itkImageProcessor::SetMetric(tMetricTypeEnum eMetric){
if(eMetric != tMetricTypeEnum::NCC &&
eMetric != tMetricTypeEnum::MI){
itkExceptionMacro(<< "Unkown metric string : " << eMetric);
return;
}
std::cout<<"Setting Metric: "<<eMetric<<std::endl;
m_r23MetaInfo->SetMetricType(eMetric);
}
void itkImageProcessor::SetOptimizer(std::string optimizer)
{
@ -1167,11 +1312,59 @@ void itkImageProcessor::SetMetric(std::string metric)
}
void itkImageProcessor::SetMaxNumberOfIterations(int iNum)
{
m_r23MetaInfo->SetMaxIterations(iNum);
void itkImageProcessor::SetDegreeOfFreedom(tDegreeOfFreedomEnum dof){
m_DRTGeometryMetaInfo->SetDegreeOfFreedom(dof);
m_r23MetaInfo->SetDegreeOfFreedom(dof);
}
void itkImageProcessor::SetHandleRotationImportOffset(eHandlingRotImpTransform hrio){
m_DRTGeometryMetaInfo->SetHandleRotationImpOffset(hrio);
}
void itkImageProcessor::SetPowellOptimParameters(
double dStepT,
double dValTol,
double dStepL,
int iMaxLinI,
int iMaxIter){
m_PowellOptimizerMetaInfo->SetStepTolerance(dStepT);
m_PowellOptimizerMetaInfo->SetValueTolerance(dValTol);
m_PowellOptimizerMetaInfo->SetStepLength(dStepL);
m_PowellOptimizerMetaInfo->SetMaximumLineInteration(iMaxLinI);
m_PowellOptimizerMetaInfo->SetMaxIterations(iMaxIter);
}
void itkImageProcessor::SetAmoebaOptimParameters(
double dParConvTol,
double dFuntConvTol,
double dSimplex,
int iMaxIter){
m_AmoebaOptimizerMetaInfo->SetParametersConvergenceTolerance(dParConvTol);
m_AmoebaOptimizerMetaInfo->SetFunctionConvergenceTolerance(dFuntConvTol);
m_AmoebaOptimizerMetaInfo->SetSimplexDelta(dSimplex);
m_AmoebaOptimizerMetaInfo->SetMaxIterations(iMaxIter);
}
void itkImageProcessor::SetMIMetricParameters(double dMaxT,int iNhb){
m_MIMetricMetaInfo->SetMaxTranslation(dMaxT);
m_MIMetricMetaInfo->SetNumberOfHistogramBins(iNhb);
}
void itkImageProcessor::SetNCCMetricParameters(double dMaxT,bool bSm){
m_NCCMetricMetaInfo->SetMaxTranslation(dMaxT);
m_NCCMetricMetaInfo->SetSubtractMean(bSm);
}
void itkImageProcessor::InitializeProjector()
{
@ -1179,7 +1372,7 @@ void itkImageProcessor::InitializeProjector()
m_DRTImage2MetaInfo == NULL ||
m_DRTGeometryMetaInfo == NULL ||
m_TransformMetaInfo == NULL ){
return;
//TODO
return;
}
@ -1190,12 +1383,6 @@ void itkImageProcessor::InitializeProjector()
ImageType3D::PointType ZeroPoint;
ZeroPoint.Fill(0.);
// InternalImageType::DirectionType IECtoLPS_Directions;
// IECtoLPS_Directions =
// m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
TransformType::Pointer CurrTransform;
CurrTransform= CalculateInternalTransformV3(
@ -1220,7 +1407,6 @@ void itkImageProcessor::InitializeProjector()
transform1->SetCenter(
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() );
// 2D Image 1
interpolator1->SetProjectionAngle(
dtr * m_DRTImage1MetaInfo->GetProjectionAngleLPS() );
@ -1229,7 +1415,7 @@ void itkImageProcessor::InitializeProjector()
interpolator1->SetThreshold(
m_DRTGeometryMetaInfo->GetIntensityThreshold() );
interpolator1->SetPanelOffset(
m_DRTGeometryMetaInfo->GetPanel1Offset());
m_DRTImage1MetaInfo->GetPanelOffsetPGeo());
interpolator1->SetTransform(transform1);
interpolator1->Initialize();
@ -1253,10 +1439,8 @@ void itkImageProcessor::InitializeProjector()
CurrTransform->GetAngleZ()
);
transform2->SetCenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() );
//transform2->Print(std::cout);
// 2D Image 2
interpolator2->SetProjectionAngle(
@ -1266,7 +1450,7 @@ void itkImageProcessor::InitializeProjector()
interpolator2->SetThreshold(
m_DRTGeometryMetaInfo->GetIntensityThreshold() );
interpolator2->SetPanelOffset(
m_DRTGeometryMetaInfo->GetPanel2Offset());
m_DRTImage2MetaInfo->GetPanelOffsetPGeo());
interpolator2->SetTransform(transform2);
interpolator2->Initialize();
@ -1274,7 +1458,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);
@ -1293,7 +1481,8 @@ void itkImageProcessor::InitializeProjector()
resampleFilter2 = ResampleFilterType::New();
resampleFilter2->SetInput(
image3DIn);
m_VolumeSourceDupli->GetOutput());
resampleFilter2->SetDefaultPixelValue(0);
resampleFilter2->SetNumberOfWorkUnits(iNWUnits);
@ -1491,10 +1680,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;
@ -1505,17 +1694,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(
@ -1541,19 +1730,17 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
m_DRTGeometryMetaInfo->GetDRT1Spacing()
);
ImageType3D::PointType PanelOffsetIEC1;
PanelOffsetIEC1[0] = -m_DRTGeometryMetaInfo->GetPanel1Offset();
PanelOffsetIEC1[1] = 0.0;
PanelOffsetIEC1[2] = 0.0;
m_DRTImage1MetaInfo->SetPanelOffsetPGeo(
this->CalcPanelOffsetPGeo(m_CTMetaInfo->GetPatientOrientation(),
-m_DRTGeometryMetaInfo->GetPanel1Offset())
);
// This HAS TO be calculated from the nominal isocenter.
m_DRTImage1MetaInfo->SetOriginLPS(
CalcDRTImageOrigin(
m_DRTImage1MetaInfo->GetOrigin(),
m_DRTImage1MetaInfo->GetCOV(),
NominalIsocenter_wZcorrectionLPS,//m_DRTImage1MetaInfo->GetProjectionOriginLPS(),
NominalIsocenter_wZcorrectionLPS,
this->CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() ),
@ -1641,19 +1828,15 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
m_DRTGeometryMetaInfo->GetDRT2Spacing()
);
ImageType3D::PointType PanelOffsetIEC2;
PanelOffsetIEC2[0] = -m_DRTGeometryMetaInfo->GetPanel2Offset();
PanelOffsetIEC2[1] = 0.0;
PanelOffsetIEC2[2] = 0.0;
m_DRTImage2MetaInfo->SetPanelOffsetPGeo(
this->CalcPanelOffsetPGeo(m_CTMetaInfo->GetPatientOrientation(),
-m_DRTGeometryMetaInfo->GetPanel2Offset()));
m_DRTImage2MetaInfo->SetOriginLPS(
CalcDRTImageOrigin(
m_DRTImage2MetaInfo->GetOrigin(),
m_DRTImage2MetaInfo->GetCOV(),
NominalIsocenter_wZcorrectionLPS,//m_DRTImage2MetaInfo->GetProjectionOriginLPS(),
NominalIsocenter_wZcorrectionLPS,
this->CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) ,
@ -1661,19 +1844,6 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
));
m_DRTImage2MetaInfo->SetOriginLPS(
CalcDRTImageOrigin(
m_DRTImage2MetaInfo->GetOrigin(),
m_DRTImage2MetaInfo->GetCOV(),
NominalIsocenter_wZcorrectionLPS,//m_DRTImage2MetaInfo->GetProjectionOriginLPS(),
this->CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) ,
Std_DRT2LPS
)
);
m_DRTImage2MetaInfo->SetImageDirectionsLPS(
CalcDRTImageDirections(
this->CalcProjectionAngleLPS(
@ -1818,6 +1988,8 @@ void itkImageProcessor::GetProjectionImages(){
CurrTransform->GetAngleZ()
);
transform2 ->SetCenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero());
@ -2381,13 +2553,7 @@ itkImageProcessor::GetUserIsocentricTransform(){
TransformType::Pointer
itkImageProcessor::GetCompleteIsocentricTransform
(/*ImageType3D::PointType TransformCenter*/){
TransformType::Pointer
regTransform = TransformType::New();
regTransform->SetComputeZYX(true);
regTransform->SetIdentity();
itkImageProcessor::GetCompleteIsocentricTransform(){
if(m_TransformMetaInfo == nullptr ||
m_CTMetaInfo == nullptr){
@ -2479,23 +2645,28 @@ void itkImageProcessor::SetRegionFixed2(
void itkImageProcessor::ResetROIRegions(){
auto region1temp = m_PASourceDupli->GetOutput()->GetBufferedRegion();
auto region2temp = m_LATSourceDupli->GetOutput()->GetBufferedRegion();
// auto region1temp = m_PASourceDupli->GetOutput()->GetBufferedRegion();
// auto region2temp = m_LATSourceDupli->GetOutput()->GetBufferedRegion();
this->m_R23->SetroiAutoReg1(region1temp);
this->m_R23->SetroiAutoReg2(region2temp);
// this->m_R23->SetroiAutoReg1(region1temp);
// this->m_R23->SetroiAutoReg2(region2temp);
this->m_R23->setROISizeToZero();
}
void itkImageProcessor::InizializeRegistration(double stepLength,
double maxTranslation,
eDegreeOfFreedomType dof){
void itkImageProcessor::InizializeRegistration(){
m_r23MetaInfo->SetDegreeOfFreedom(dof);
m_R23->SetPA(m_PASourceDupli->GetOutput());
m_R23->SetLAT(m_LATSourceDupli->GetOutput());
m_R23->SetVolume(m_VolumeSourceDupli->GetOutput());
m_R23->Setr23Meta(m_r23MetaInfo);
m_R23->SetPowellMeta(m_PowellOptimizerMetaInfo);
m_R23->SetAmoebaMeta(m_AmoebaOptimizerMetaInfo);
m_R23->SetMIMeta(m_MIMetricMetaInfo);
m_R23->SetNCCMeta(m_NCCMetricMetaInfo);
m_R23->SetInternalTransf1(m_InternalTransf1);
m_R23->SetInternalTransf2(m_InternalTransf2);
m_R23->Setfilter1(filter1);
@ -2504,8 +2675,7 @@ void itkImageProcessor::InizializeRegistration(double stepLength,
m_R23->Setinterpolator2(interpolator2);
m_R23->SetTransformMetaInfo(m_TransformMetaInfo);
m_R23->InitializeRegistration(
stepLength,maxTranslation,dof);
m_R23->InitializeRegistration();
}

67
reg23Topograms/itkDTRrecon/itkImageProcessor.h
View File

@ -47,6 +47,7 @@ gfattori 08.11.2021
#include "itkImageProcessorHelpers.h"
#include "itkReg23.h"
#include "itkReg23MetaInformation.h"
namespace itk
{
@ -70,6 +71,12 @@ public:
itkTypeMacro(itkImageProcessor, Object);
itkSetMacro(ResamplepCT, bool);
itkGetMacro(ResamplepCT, bool);
itkSetMacro(SamplingLNG, double);
itkGetMacro(SamplingLNG, double);
CommandIterationUpdate::Pointer
optimizerObserver;
@ -86,19 +93,37 @@ public:
void SetRegionFixed1(int,int,int,int);
void SetRegionFixed2(int,int,int,int);
void ResetROIRegions();
void InizializeRegistration(double stepLength,
double maxTranslation,
eDegreeOfFreedomType dof);
void InizializeRegistration();
/** Maximum number of iterations for the optimizer */
void SetMaxNumberOfIterations(int);
itkGetMacro(R23, itkReg23::Pointer);
/** Set Optimizer type */
void SetOptimizer(std::string);
void SetOptimizer(tOptimizerTypeEnum);
/** Set Metric type */
void SetMetric(std::string);
void SetMetric(tMetricTypeEnum);
/** Set DOF */
void SetDegreeOfFreedom(tDegreeOfFreedomEnum);
/** Set Handle Of Rotation Import Offset */
void SetHandleRotationImportOffset(tHandlingRotImpTransform);
/** Set PowellOptimizer */
void SetPowellOptimParameters(double,double,double,int,int);
/** Set AmoebaOptimizer */
void SetAmoebaOptimParameters(double,double,double,int);
/** Set MI */
void SetMIMetricParameters(double,int);
/** Set NCC */
void SetNCCMetricParameters(double,bool);
/** Set number of logic CPU to be made available to interpolators*/
void SetNumberOfWorkingUnits(int iN);
@ -133,6 +158,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();
@ -238,6 +267,11 @@ protected:
virtual ~itkImageProcessor();
void PrintSelf(std::ostream& os, itk::Indent indent) const;
bool
m_ResamplepCT;
double
m_SamplingLNG;
private:
itkImageProcessor(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
@ -267,6 +301,11 @@ private:
using ImageIOType = itk::GDCMImageIO;
using DuplicatorType = itk::ImageDuplicator<InternalImageType>;
/** Image input resampling Types */
using ResampleInputFilterType = itk::ResampleImageFilter<ImageType3D, ImageType3D>;
using LinearInterpolatorType = itk::LinearInterpolateImageFunction<ImageType3D, double>;
using SizeValueType = ImageType3D::SizeType::SizeValueType;
/** widely used filters: flip, cast, rescale, ... */
using Input2DRescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, InternalImageType>;
using CastFilterType3D = itk::CastImageFilter<ImageType3D, InternalImageType>;
@ -292,7 +331,6 @@ private:
toVTKLocalizer2;
InternalImageType::Pointer
image3DIn,
imageDRT1In,
imageDRT2In;
@ -303,9 +341,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. */
@ -333,6 +368,10 @@ private:
tPatOrientation pOrient,
double pAngleIEC);
double
CalcPanelOffsetPGeo(
tPatOrientation pOrient,
double dPOffset);
/** Apply transform to CT volume.
* This moves the volume based on RT Plan info
* if those are available. */
@ -387,6 +426,18 @@ private:
R23MetaInformation::Pointer
m_r23MetaInfo;
PowellOptimizerMetaInformation::Pointer
m_PowellOptimizerMetaInfo;
AmoebaOptimizerMetaInformation::Pointer
m_AmoebaOptimizerMetaInfo;
MIMetricMetaInformation::Pointer
m_MIMetricMetaInfo;
NCCMetricMetaInformation::Pointer
m_NCCMetricMetaInfo;
InternalTransformMetaInformation::Pointer
m_InternalTransf1,
m_InternalTransf2;

15
reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp
View File

@ -100,6 +100,7 @@ CalculateInternalTransformV3(
const char* queryStationName(const char* pcFName){
static std::string buffer;
buffer.clear();
/* Check if we can open the file */
gdcm::Reader R;
@ -112,10 +113,8 @@ const char* queryStationName(const char* pcFName){
const gdcm::File &file = R.GetFile();
const gdcm::DataSet &ds = file.GetDataSet();
std::string s (gGetStringValueFromTag(gdcm::Tag(0x0008,0x1010), ds));
std::copy(s.rbegin(), s.rend(), std::back_inserter(buffer));
std::reverse(buffer.begin(),buffer.end());
std::string s(gGetStringValueFromTag(gdcm::Tag(0x0008,0x1010), ds));
std::copy(s.begin(), s.end(), buffer.begin());
return buffer.c_str();
@ -134,10 +133,12 @@ int query2DimageReconstructionDiameter(const char *pcFName){
const gdcm::File &file = R.GetFile();
const gdcm::DataSet &ds = file.GetDataSet();
char* sTmpString = new char [255];
strcpy( sTmpString,gGetStringValueFromTag(gdcm::Tag(0x0018,0x1100), ds));
std::string buffer;
return std::atoi(sTmpString);
std::string s(gGetStringValueFromTag(gdcm::Tag(0x0018,0x1100), ds));
std::copy(s.begin(), s.end(), buffer.begin());
return std::stoi(buffer);
}

2
reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h
View File

@ -49,7 +49,7 @@ int query2DimageReconstructionDiameter(const char*);
const char* queryStationName(const char*);
}
};
#endif

84
reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h
View File

@ -28,9 +28,25 @@ public slots:
void onAbortProcessCommand(){
bAbortProcessCommand = true;
};
void onIteration(double dI,double dX,double dY,double dZ){
emit
sendRegistrationProgress(dI,dX,dY,dZ);
void onIteration(itk::tOptimizerTypeEnum eType,double dI,double dX,double dY,double dZ){
switch (eType) {
case itk::eOptimizerType::POWELL:
emit
sendRegistrationProgress(0,dI,dX,dY,dZ);
break;
case itk::eOptimizerType::AMOEBA:
emit
sendRegistrationProgress(1,dI,dX,dY,dZ);
break;
default:
break;
}
}
private:
@ -42,7 +58,7 @@ private:
signals:
void sendRegistrationProgress(double,double,double,double);
void sendRegistrationProgress(int,double,double,double,double);
};
@ -57,6 +73,11 @@ class CommandIterationUpdate : public itk::Command {
public:
QObjectIterationUpdate *objIterUpdate;
tOptimizerTypeEnum
eOptimType;
int
iAmoebaIterations;
using Self = CommandIterationUpdate;
using Superclass = itk::Command;
@ -85,7 +106,12 @@ public:
using OptimizerPointer = const OptimizerType*;
using AmoebaOptimizerType = itk::AmoebaOptimizer;
using AmoebaOptimizerPointer = const OptimizerType*;
using AmoebaOptimizerPointer = const AmoebaOptimizerType*;
void setOptimizerTypeOfIterationUpdate(tOptimizerTypeEnum eOptType){
eOptimType = eOptType;
iAmoebaIterations = 0;
}
void
Execute(itk::Object* caller, const itk::EventObject& event) override
@ -103,34 +129,40 @@ public:
}
// std::cout << "Progress: " << this->m_Process->GetAbortGenerateData() << std::endl;
auto optimizer = dynamic_cast<OptimizerPointer>(object);
OptimizerPointer optPow;
AmoebaOptimizerPointer optAm;
if (typeid(event) == typeid(itk::IterationEvent)) {
switch (eOptimType) {
case eOptimizerType::POWELL:
optPow = dynamic_cast<OptimizerPointer>(object);
if (typeid(event) == typeid(itk::IterationEvent)) {
//Feedback from the optimizer executed at the end of every itteration
// currently just print the result into the cout. We might add
// functionality to register and emit signals to update the UI.
objIterUpdate->onIteration(
eOptimType,
optPow->GetCurrentIteration()+1,
optPow->GetCurrentPosition()[0],
optPow->GetCurrentPosition()[1],
optPow->GetCurrentPosition()[2]
);
}
std::cout << "Iteration: " << optimizer->GetCurrentIteration() << std::endl;
auto oldprecision = std::cout.precision();
std::cout.precision(std::numeric_limits<double>::digits10 + 2);
std::cout << "Similarity: " << optimizer->GetValue() << std::endl;
std::cout.precision(oldprecision);
std::cout << "Position: " << optimizer->GetCurrentPosition() << std::endl;
break;
case eOptimizerType::AMOEBA:
optAm = dynamic_cast<AmoebaOptimizerPointer>(object);
iAmoebaIterations ++;
// objIterUpdate->onIteration(
// optimizer->GetCurrentIteration()+1,
// optimizer->GetCurrentPosition()[0],
// optimizer->GetCurrentPosition()[2],
// -optimizer->GetCurrentPosition()[1]
// );
objIterUpdate->onIteration(
optimizer->GetCurrentIteration()+1,
optimizer->GetCurrentPosition()[0],
optimizer->GetCurrentPosition()[1],
optimizer->GetCurrentPosition()[2]
eOptimType,
iAmoebaIterations, //optAm->GetCachedValue(),
optAm->GetCachedCurrentPosition()[0],
optAm->GetCachedCurrentPosition()[1],
optAm->GetCachedCurrentPosition()[2]
);
break;
default:
break;
}
return;
}
};

102
reg23Topograms/itkDTRrecon/itkReg23.cpp
View File

@ -13,6 +13,7 @@ itkReg23::itkReg23()
NCCmetric = MetricType::New();
MImetric = MIMetricType::New();
// Optimizers
PowellOptimizer = PowellOptimizerType::New();
AmoebaOptimizer = AmoebaOptimizerType::New();
@ -46,6 +47,8 @@ itkReg23::itkReg23()
// m_UseFullROI = true; // if the full image ROI shall be used
m_UseDumptoFile = false; //If each (!) optimzer result shall be dumpped into a file.
this->setROISizeToZero();
}
@ -54,15 +57,22 @@ itkReg23::~itkReg23()
}
void itkReg23::setROISizeToZero(){
ImageType3D::SizeType zeroSize;
zeroSize.Fill(0);
m_roiAutoReg1.SetSize(zeroSize);
m_roiAutoReg2.SetSize(zeroSize);
}
void itkReg23::PrintSelf(std::ostream& os, Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
void itkReg23::InitializeRegistration(
double stepLength,
double maxTranslation,
eDegreeOfFreedomType dof)
void itkReg23::InitializeRegistration()
{
std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
@ -114,9 +124,19 @@ void itkReg23::InitializeRegistration(
itkExceptionMacro(<< "m_interpolator2 data not present");
}
AmoebaOptimizerType::ParametersType simplexDelta(3);
ExhaustiveOptimizerType::StepsType steps(3);
const int numberOfSteps = 25; //In each direction. Total number of steps is ((2*numberOfSteps+1))^3. For 25 -> 132651.
// std::cout << "Print r23" <<
// m_r23Meta->GetMetricType() << " " << m_r23Meta->GetDegreeOfFreedom() << std::endl;
// std::cout << "Print Powell" <<
// m_PowellMeta->GetStepLength() << " " << m_PowellMeta->GetStepTolerance() << std::endl;
// std::cout << "Print Amoeba" <<
// m_AmoebaMeta->GetFunctionConvergenceTolerance() << " " << m_AmoebaMeta->GetSimplexDelta() << std::endl;
// std::cout << "Print NCC" <<
// m_NCCMeta->GetSubtractMean() << " " << m_NCCMeta->GetMaxTranslation() << std::endl;
// std::cout << "Print MI" <<
// m_MIMeta->GetMaxTranslation() << std::endl;
// ExhaustiveOptimizerType::StepsType steps(3);
// const int numberOfSteps = 25; //In each direction. Total number of steps is ((2*numberOfSteps+1))^3. For 25 -> 132651.
//const double stepLength = 0.1;
switch (m_r23Meta->GetOptimizerType()) {
@ -125,13 +145,14 @@ void itkReg23::InitializeRegistration(
std::cout<< "Using POWELL Optimizer" <<std::endl;
PowellOptimizer->SetMaximize(false); // for NCC and MI
PowellOptimizer->SetMaximumIteration(m_r23Meta->GetMaxIterations());
PowellOptimizer->SetMaximumLineIteration(4); // for Powell's method
PowellOptimizer->SetStepLength(stepLength);
PowellOptimizer->SetStepTolerance(0.01);
PowellOptimizer->SetValueTolerance(0.000002);
PowellOptimizer->SetMaximumIteration(m_PowellMeta->GetMaxIterations());
PowellOptimizer->SetMaximumLineIteration(m_PowellMeta->GetMaximumLineInteration()); // for Powell's method
PowellOptimizer->SetStepLength(m_PowellMeta->GetStepLength());
PowellOptimizer->SetStepTolerance(m_PowellMeta->GetStepTolerance());
PowellOptimizer->SetValueTolerance(m_PowellMeta->GetValueTolerance());
PowellOptimizer->RemoveAllObservers();
PowellOptimizer->AddObserver(itk::IterationEvent(), m_OptimizerObserver);
m_OptimizerObserver->setOptimizerTypeOfIterationUpdate(tOptimizerTypeEnum::POWELL);
registration->SetOptimizer(PowellOptimizer);
@ -141,21 +162,33 @@ void itkReg23::InitializeRegistration(
std::cout<< "Using AMOEBA Optimizer" <<std::endl;
AmoebaOptimizer->SetMinimize(true);
AmoebaOptimizer->SetMaximumNumberOfIterations(m_r23Meta->GetMaxIterations());
AmoebaOptimizer->SetParametersConvergenceTolerance(0.1);
AmoebaOptimizer->SetFunctionConvergenceTolerance(0.000002);
AmoebaOptimizer->SetMaximumNumberOfIterations(m_AmoebaMeta->GetMaxIterations());
AmoebaOptimizer->SetParametersConvergenceTolerance(m_AmoebaMeta->GetParametersConvergenceTolerance());
AmoebaOptimizer->SetFunctionConvergenceTolerance(m_AmoebaMeta->GetFunctionConvergenceTolerance());
AmoebaOptimizer->SetAutomaticInitialSimplex(false);
//Initial size of the simplex (otherwise it is a very small number and optimizer stops immeditaly)
// 2 mm / 2 degrees seems to be a good value which performs well.
if (GetNumberOfDOF(dof) == 0) {
if (GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()) == 0) {
itkExceptionMacro(<< "Unkown or unsupported degree of freedom");
return;
}else if (GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()) == 3){
AmoebaOptimizerType::ParametersType simplexDelta3dof(3);
for (int i = 0; i < GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()); i++) {
simplexDelta3dof[i] = m_AmoebaMeta->GetSimplexDelta();
}
AmoebaOptimizer->SetInitialSimplexDelta(simplexDelta3dof);
}else if (GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()) == 6){
AmoebaOptimizerType::ParametersType simplexDelta6dof(6);
for (int i = 0; i < GetNumberOfDOF(m_r23Meta->GetDegreeOfFreedom()); i++) {
simplexDelta6dof[i] = m_AmoebaMeta->GetSimplexDelta();
}
AmoebaOptimizer->SetInitialSimplexDelta(simplexDelta6dof);
}
for (int i = 0; i < GetNumberOfDOF(dof); i++) {
simplexDelta[i] = 2.0;
}
AmoebaOptimizer->SetInitialSimplexDelta(simplexDelta);
AmoebaOptimizer->RemoveAllObservers();
AmoebaOptimizer->AddObserver(itk::IterationEvent(), m_OptimizerObserver);
m_OptimizerObserver->setOptimizerTypeOfIterationUpdate(tOptimizerTypeEnum::AMOEBA);
registration->SetOptimizer(AmoebaOptimizer);
@ -164,13 +197,13 @@ void itkReg23::InitializeRegistration(
case tOptimizerTypeEnum::EXHAUSTIVE:
std::cout<< "Using Extensive Optimizer" <<std::endl;
steps[0] = numberOfSteps;
steps[1] = numberOfSteps;
steps[2] = numberOfSteps;
ExhaustiveOptimizer->SetNumberOfSteps(steps);
ExhaustiveOptimizer->SetStepLength(stepLength);
// steps[0] = numberOfSteps;
// steps[1] = numberOfSteps;
// steps[2] = numberOfSteps;
// ExhaustiveOptimizer->SetNumberOfSteps(steps);
// ExhaustiveOptimizer->SetStepLength(stepLength);
registration->SetOptimizer(ExhaustiveOptimizer);
// registration->SetOptimizer(ExhaustiveOptimizer);
break;
@ -185,8 +218,9 @@ void itkReg23::InitializeRegistration(
std::cout<< "Using MI Metric" <<std::endl;
registration->SetMetric(MImetric);
MImetric->SetNumberOfHistogramBins(m_MIMeta->GetNumberOfHistogramBins());
MImetric->ComputeGradientOff();
MImetric->SetMaxTranslation(maxTranslation);
MImetric->SetMaxTranslation(m_MIMeta->GetMaxTranslation());
break;
@ -195,7 +229,7 @@ void itkReg23::InitializeRegistration(
registration->SetMetric(NCCmetric);
NCCmetric->ComputeGradientOff();
NCCmetric->SetSubtractMean(true);
NCCmetric->SetMaxTranslation(maxTranslation);
NCCmetric->SetMaxTranslation(m_NCCMeta->GetMaxTranslation());
break;
@ -208,6 +242,18 @@ void itkReg23::InitializeRegistration(
registration->SetFixedImage2(m_LAT);
registration->SetMovingImage(m_Volume);
if( this->GetroiAutoReg1().GetNumberOfPixels() == 0 ){
auto defaultRegion = m_PA->GetBufferedRegion();
this->SetroiAutoReg1(defaultRegion);
std::cout<< "Image1 region not defined. Using the whole image buffered region" <<std::endl;
}
if( this->GetroiAutoReg2().GetNumberOfPixels() == 0 ){
auto defaultRegion = m_LAT->GetBufferedRegion();
this->SetroiAutoReg2(defaultRegion);
std::cout<< "Image2 region not defined. Using the whole image buffered region" <<std::endl;
}
registration->SetFixedImageRegion1(m_roiAutoReg1);
registration->SetFixedImageRegion2(m_roiAutoReg2);

33
reg23Topograms/itkDTRrecon/itkReg23.h
View File

@ -18,6 +18,8 @@
#include "itkImageProcessorHelpers.h"
#include "itkQtIterationUpdate.h"
#include "itkReg23MetaInformation.h"
namespace itk{
@ -50,6 +52,19 @@ public:
itkSetMacro(r23Meta, R23MetaInformation::Pointer);
itkGetMacro(r23Meta, R23MetaInformation::Pointer);
itkSetMacro(PowellMeta, PowellOptimizerMetaInformation::Pointer);
itkGetMacro(PowellMeta, PowellOptimizerMetaInformation::Pointer);
itkSetMacro(AmoebaMeta, AmoebaOptimizerMetaInformation::Pointer);
itkGetMacro(AmoebaMeta, AmoebaOptimizerMetaInformation::Pointer);
itkSetMacro(MIMeta, MIMetricMetaInformation::Pointer);
itkGetMacro(MIMeta, MIMetricMetaInformation::Pointer);
itkSetMacro(NCCMeta, NCCMetricMetaInformation::Pointer);
itkGetMacro(NCCMeta, NCCMetricMetaInformation::Pointer);
itkSetMacro(Volume, InternalImageType::Pointer);
itkGetMacro(Volume, InternalImageType::Pointer);
@ -92,9 +107,7 @@ public:
/** Auto Reg23 methods */
/** Initialize the registration pipeline*/
void InitializeRegistration(double stepLength,
double maxTranslation,
eDegreeOfFreedomType dof);
void InitializeRegistration();
/** Start the registration process*/
int StartRegistration(std::string extraInfo);
/** Get the current cost function value from the optimizer*/
@ -108,6 +121,8 @@ public:
/** Auto Reg23 methods END */
void setROISizeToZero();
protected:
itkReg23();
virtual ~itkReg23();
@ -149,6 +164,18 @@ private:
R23MetaInformation::Pointer
m_r23Meta;
PowellOptimizerMetaInformation::Pointer
m_PowellMeta;
AmoebaOptimizerMetaInformation::Pointer
m_AmoebaMeta;
MIMetricMetaInformation::Pointer
m_MIMeta;
NCCMetricMetaInformation::Pointer
m_NCCMeta;
InternalImageType::Pointer
m_Volume,
m_PA,

104
reg23Topograms/itkDTRrecon/itkReg23MetaInformation.cpp Normal file
View File

@ -0,0 +1,104 @@
#include "itkReg23MetaInformation.h"
namespace itk {
PowellOptimizerMetaInformation::
PowellOptimizerMetaInformation(){
this->m_MaximumLineInteration = 4;
this->m_StepLength = 2.0;
this->m_StepTolerance = 0.01;
this->m_ValueTolerance = 0.000002;
this->m_MaxIterations = 6;
}
void
PowellOptimizerMetaInformation
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
PowellOptimizerMetaInformation
::~PowellOptimizerMetaInformation ()
{
}
AmoebaOptimizerMetaInformation::
AmoebaOptimizerMetaInformation(){
this->m_ParametersConvergenceTolerance = 0.1;
this->m_FunctionConvergenceTolerance = 0.000002;
this->m_SimplexDelta = 2.0;
this->m_MaxIterations = 100;
}
void
AmoebaOptimizerMetaInformation
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
AmoebaOptimizerMetaInformation
::~AmoebaOptimizerMetaInformation ()
{
}
MIMetricMetaInformation::
MIMetricMetaInformation(){
this->m_MaxTranslation = 100.0;
this->m_NumberOfHistogramBins = 50;
}
void
MIMetricMetaInformation
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
MIMetricMetaInformation
::~MIMetricMetaInformation ()
{
}
NCCMetricMetaInformation::
NCCMetricMetaInformation(){
this->m_MaxTranslation = 100.0;
this->m_SubtractMean = true;
}
void
NCCMetricMetaInformation
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
NCCMetricMetaInformation
::~NCCMetricMetaInformation ()
{
}
}

232
reg23Topograms/itkDTRrecon/itkReg23MetaInformation.h Normal file
View File

@ -0,0 +1,232 @@
#ifndef ITKREG23METAINFORMATION_H
#define ITKREG23METAINFORMATION_H
#include "itkObject.h"
#include "itkObjectFactory.h"
#include "itkSmartPointer.h"
#include "itkMacro.h"
namespace itk {
class PowellOptimizerMetaInformation :
public itk::Object{
public:
/** standard typedefs **/
typedef PowellOptimizerMetaInformation Self;
typedef itk::Object Superclass;
typedef itk::SmartPointer<Self> Pointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(PowellOptimizerMetaInformation, itk::Object);
/** object information streaming **/
void PrintSelf(std::ostream& os, itk::Indent indent) const;
itkSetMacro(StepTolerance,double);
itkGetMacro(StepTolerance,double);
itkSetMacro(ValueTolerance,double);
itkGetMacro(ValueTolerance,double);
itkSetMacro(StepLength,double);
itkGetMacro(StepLength,double);
itkSetMacro(MaximumLineInteration,int);
itkGetMacro(MaximumLineInteration,int);
itkSetEnumMacro(MaxIterations, int);
itkGetEnumMacro(MaxIterations, int);
protected:
double m_StepTolerance;
double m_ValueTolerance;
double m_StepLength;
int m_MaximumLineInteration;
int
m_MaxIterations;
/** Default Constructor **/
PowellOptimizerMetaInformation ();
/** Default Destructor **/
virtual ~PowellOptimizerMetaInformation ();
private:
/** purposely not implemented **/
PowellOptimizerMetaInformation (const Self&);
/** purposely not implemented **/
void operator=(const Self&);
};
class AmoebaOptimizerMetaInformation :
public itk::Object{
public:
/** standard typedefs **/
typedef AmoebaOptimizerMetaInformation Self;
typedef itk::Object Superclass;
typedef itk::SmartPointer<Self> Pointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(AmoebaOptimizerMetaInformation, itk::Object);
/** object information streaming **/
void PrintSelf(std::ostream& os, itk::Indent indent) const;
itkSetMacro(ParametersConvergenceTolerance,double);
itkGetMacro(ParametersConvergenceTolerance,double);
itkSetMacro(FunctionConvergenceTolerance,double);
itkGetMacro(FunctionConvergenceTolerance,double);
itkSetMacro(SimplexDelta,double);
itkGetMacro(SimplexDelta,double);
itkSetEnumMacro(MaxIterations, int);
itkGetEnumMacro(MaxIterations, int);
protected:
double m_ParametersConvergenceTolerance;
double m_FunctionConvergenceTolerance;
double m_SimplexDelta;
int
m_MaxIterations;
/** Default Constructor **/
AmoebaOptimizerMetaInformation ();
/** Default Destructor **/
virtual ~AmoebaOptimizerMetaInformation ();
private:
/** purposely not implemented **/
AmoebaOptimizerMetaInformation (const Self&);
/** purposely not implemented **/
void operator=(const Self&);
};
class MIMetricMetaInformation :
public itk::Object{
public:
/** standard typedefs **/
typedef MIMetricMetaInformation Self;
typedef itk::Object Superclass;
typedef itk::SmartPointer<Self> Pointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(MIMetricMetaInformation, itk::Object);
/** object information streaming **/
void PrintSelf(std::ostream& os, itk::Indent indent) const;
itkSetMacro(MaxTranslation,double);
itkGetMacro(MaxTranslation,double);
itkSetMacro(NumberOfHistogramBins,int);
itkGetMacro(NumberOfHistogramBins,int);
protected:
double m_MaxTranslation;
int m_NumberOfHistogramBins;
/** Default Constructor **/
MIMetricMetaInformation ();
/** Default Destructor **/
virtual ~MIMetricMetaInformation ();
private:
/** purposely not implemented **/
MIMetricMetaInformation (const Self&);
/** purposely not implemented **/
void operator=(const Self&);
};
class NCCMetricMetaInformation :
public itk::Object{
public:
/** standard typedefs **/
typedef NCCMetricMetaInformation Self;
typedef itk::Object Superclass;
typedef itk::SmartPointer<Self> Pointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(NCCMetricMetaInformation, itk::Object);
/** object information streaming **/
void PrintSelf(std::ostream& os, itk::Indent indent) const;
itkSetMacro(MaxTranslation,double);
itkGetMacro(MaxTranslation,double);
itkSetMacro(SubtractMean,bool);
itkGetMacro(SubtractMean,bool);
protected:
double m_MaxTranslation;
bool m_SubtractMean;
/** Default Constructor **/
NCCMetricMetaInformation ();
/** Default Destructor **/
virtual ~NCCMetricMetaInformation ();
private:
/** purposely not implemented **/
NCCMetricMetaInformation (const Self&);
/** purposely not implemented **/
void operator=(const Self&);
};
}
#endif

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.;
PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);
// Get ths input pointers
InputImageConstPointer inputPtr = this->GetInputImage();
//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

34
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;
@ -110,7 +111,7 @@ int vtkContourTopogramProjectionFilter::RequestData(
// If no points, then nothing to do.
if (points == nullptr)
{
std::cout << "Cannot Project; no input points" << std::endl;
// std::cout << "Cannot Project; no input points" << std::endl;
vtkDebugMacro("Cannot Project; no input points");
return 1;
}
@ -118,7 +119,7 @@ int vtkContourTopogramProjectionFilter::RequestData(
// If reference transform, then nothing to do.
if (m_RefTransform == nullptr)
{
std::cout << "Cannot Project; no input reference projection transform" << std::endl;
// std::cout << "Cannot Project; no input reference projection transform" << std::endl;
vtkDebugMacro("Cannot Project; no input reference projection transform");
return 1;
}
@ -126,7 +127,7 @@ int vtkContourTopogramProjectionFilter::RequestData(
// If transform, then nothing to do.
if (m_Transform == nullptr)
{
std::cout << "Cannot Project; no input projection transform" << std::endl;
// std::cout << "Cannot Project; no input projection transform" << std::endl;
vtkDebugMacro("Cannot Project; no input projection transform");
return 1;
}
@ -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.");
}
@ -167,6 +168,26 @@ int vtkContourTopogramProjectionFilter::RequestData(
pp2[1] -= dProjectionLPSOffset [1];
pp2[2] -= dProjectionLPSOffset [2];
// // Workaround to account for PA-LAT HFS-FFS
// switch (this->m_ProjectionOrientation) {
// case itk::PA:
// pp2[0] += dPOffset;
// break;
// case itk::LAT:
// switch (this->m_PatientOrientation) {
// case itk::FFS:
// pp2[1] += dPOffset;
// break;
// case itk::HFS:
// pp2[1] -= dPOffset;
// break;
// default:
// break;
// }
// break;
// default:
// break;
// }
// std::cout << ii << " " << pp2[0] << " " << pp2[1] << " " << pp2[2] << std::endl;
PrjPoints->InsertPoint(ii,pp2);
@ -213,6 +234,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));

17
reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.h
View File

@ -19,6 +19,8 @@
#include "vtkMatrix4x4.h"
#include "vtkTransform.h"
#include "itkMacro.h"
#include "DRTMetaInformation.h"
class VTKFILTERSCORE_EXPORT vtkContourTopogramProjectionFilter :
public vtkPolyDataAlgorithm
@ -37,9 +39,16 @@ public:
vtkGetObjectMacro(m_Transform, vtkTransform);
vtkGetObjectMacro(m_RefTransform, vtkTransform);
itkSetEnumMacro(ProjectionOrientation, itk::tProjOrientationType);
itkGetEnumMacro(ProjectionOrientation, itk::tProjOrientationType);
itkSetEnumMacro(PatientOrientation, itk::tPatOrientation);
itkGetEnumMacro(PatientOrientation, itk::tPatOrientation);
void SetLPStoProjectionLPSGeometryOffset (const double * );
void SetSCD(const double dVal);
void SetImportOffsetLPS(const double *);
void SetPanelOffset(const double dVal);
void cleanUpFilter();
@ -51,6 +60,11 @@ protected:
vtkInformationVector* outputVector) override;
int FillInputPortInformation(int port, vtkInformation* info) override;
itk::tProjOrientationType
m_ProjectionOrientation;
itk::tPatOrientation
m_PatientOrientation;
private:
vtkContourTopogramProjectionFilter (const vtkContourTopogramProjectionFilter &) = delete;
@ -60,7 +74,8 @@ private:
double
dProjectionLPSOffset[3],
dImportOffsetLPS[3],
dSCD;
dSCD,
dPOffset;
vtkTransform *m_RefTransform;
vtkTransform *m_Transform;