diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp index 48cd6ed..40d6ff2 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp @@ -75,7 +75,6 @@ SetPatientOrientation(tPatOrientation m_orient) DRTImageMetaInformation:: DRTImageMetaInformation(){ - this->m_RequestedSize.Fill(0.); this->m_Size.Fill(0.); @@ -83,8 +82,18 @@ DRTImageMetaInformation(){ this->m_Origin.Fill(0.); + this->m_ProjectionOriginLPS.Fill(0.); + + this->m_ProjectionOriginLPSZero.Fill(0.); + this->m_ProjectionAngleLPS = 0.; + this->m_SCD = 0.; + + this->m_ImageDirectionsLPS.SetIdentity(); + + this->m_OriginLPS.Fill(0.); + } void @@ -102,6 +111,46 @@ DRTImageMetaInformation } +DRTImageMetaInformation::PointType +DRTImageMetaInformation::GetLPStoProjectionGeoLPSOffset() +{ + + PointType pOffset; + + pOffset[0] = - m_ProjectionOriginLPS [0] + m_ProjectionOriginLPSZero [0]; + pOffset[1] = - m_ProjectionOriginLPS [1] + m_ProjectionOriginLPSZero [1]; + pOffset[2] = - m_ProjectionOriginLPS [2] + m_ProjectionOriginLPSZero [2]; + + return + pOffset; +} + +DRTImageMetaInformation ::PointType +DRTImageMetaInformation ::GetCOV() +{ + + /* calculate image size in 3D Space by finding the last voxel position */ + PointType Dest; + Dest[0]=(m_Size[0]-1) * m_Spacing [0]; + Dest[1]=(m_Size[1]-1) * m_Spacing [1]; + Dest[2]=(m_Size[2]-1) * m_Spacing [2]; + + PointType LastVoxelPosition = Dest; + + LastVoxelPosition [0] += m_Origin[0]; + LastVoxelPosition [1] += m_Origin[1]; + LastVoxelPosition [2] += m_Origin[2]; + + PointType image3DCOV; + image3DCOV[0] = (m_Origin [0] + LastVoxelPosition[0])/2.; + image3DCOV[1] = (m_Origin [1] + LastVoxelPosition[1])/2.; + image3DCOV[2] = (m_Origin [2] + LastVoxelPosition[2])/2.; + + return + image3DCOV; +} + + CTVolumeImageMetaInformation:: CTVolumeImageMetaInformation(){ @@ -118,6 +167,7 @@ CTVolumeImageMetaInformation(){ this->m_LPS2IECDirections.SetIdentity(); + this->m_ImportOffset.Fill(0.); } @@ -186,6 +236,12 @@ CTVolumeImageMetaInformation::GetCOV() image3DCOVLPS; } +CTVolumeImageMetaInformation::PointType +CTVolumeImageMetaInformation::GetOriginWOffset(){ + return + this->m_OriginLPS - this->m_ImportOffset; +} + DRTProjectionGeometryImageMetaInformation:: DRTProjectionGeometryImageMetaInformation(){ @@ -198,9 +254,17 @@ DRTProjectionGeometryImageMetaInformation(){ this->m_IntensityThreshold=0.; - this->m_ProjectionOriginLPS.Fill(0.); + this->m_DRT1Size.Fill(0.); - this->m_ProjectionOriginLPSZero.Fill(0.); + this->m_DRT2Size.Fill(0.); + + this->m_DRT1Spacing.Fill(0.); + + this->m_DRT2Spacing.Fill(0.); + + this->m_IECS2IECScannerT.Fill(0.); + + this->m_IECS2IECScannerR.Fill(0.); this->m_ProjectionCenter.Fill(0.); @@ -221,4 +285,35 @@ DRTProjectionGeometryImageMetaInformation } + + +RTGeometryMetaInformation:: +RTGeometryMetaInformation(){ + + this->m_IsocenterLPS.Fill(0.); + + this->m_IsocenterIECS.Fill(0.); + + +} + +void +RTGeometryMetaInformation +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); +} + + +RTGeometryMetaInformation +::~RTGeometryMetaInformation() +{ + +} + + + + + + } diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h index b8747e4..d997ee0 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h @@ -97,8 +97,10 @@ public: typedef itk::Object Superclass; typedef itk::SmartPointer Pointer; typedef itk::Point PointType; - typedef itk::FixedArray SpacingType; - typedef itk::FixedArray SizeType; + typedef itk::Vector SpacingType; + //typedef itk::FixedArray SizeType; + typedef itk::Size<3> SizeType; + typedef itk::Matrix DirectionType; /** Method for creation through the object factory. */ itkNewMacro(Self); @@ -109,10 +111,6 @@ public: /** object information streaming **/ void PrintSelf(std::ostream& os, itk::Indent indent) const; - - itkSetMacro(RequestedSize,SizeType); - itkGetMacro(RequestedSize,SizeType); - itkSetMacro(Size,SizeType); itkGetMacro(Size,SizeType); @@ -122,23 +120,47 @@ public: itkSetMacro(Origin,PointType); itkGetMacro(Origin,PointType); + itkSetMacro(OriginLPS,PointType); + itkGetMacro(OriginLPS,PointType); + + itkSetMacro(ProjectionOriginLPS,PointType); + itkGetMacro(ProjectionOriginLPS,PointType); + + itkSetMacro(ProjectionOriginLPSZero,PointType); + itkGetMacro(ProjectionOriginLPSZero,PointType); + itkSetMacro(ProjectionAngleLPS,double); itkGetMacro(ProjectionAngleLPS,double); + itkSetMacro(SCD,double); + itkGetMacro(SCD,double); + + itkSetMacro(ImageDirectionsLPS,DirectionType); + itkGetMacro(ImageDirectionsLPS,DirectionType); + + PointType GetCOV(); + PointType GetLPStoProjectionGeoLPSOffset(); + protected: SizeType - m_RequestedSize, m_Size; SpacingType m_Spacing; PointType - m_Origin; + m_Origin, + m_OriginLPS, + m_ProjectionOriginLPS, + m_ProjectionOriginLPSZero; double - m_ProjectionAngleLPS; + m_ProjectionAngleLPS, + m_SCD; + + DirectionType + m_ImageDirectionsLPS; /** Default Constructor **/ DRTImageMetaInformation (); @@ -164,8 +186,7 @@ public: typedef itk::SmartPointer Pointer; typedef itk::Matrix DirectionType; typedef itk::Point PointType; - typedef itk::FixedArray SpacingType; - //typedef itk::FixedArray SizeType; + typedef itk::Vector SpacingType; typedef itk::Size<3> SizeType; /** Method for creation through the object factory. */ @@ -190,6 +211,9 @@ public: itkSetMacro(OriginLPS,PointType); itkGetMacro(OriginLPS,PointType); + itkSetMacro(ImportOffset,PointType); + itkGetMacro(ImportOffset,PointType); + itkSetMacro(ImageDirections,DirectionType); itkGetMacro(ImageDirections,DirectionType); @@ -197,6 +221,7 @@ public: itkGetMacro(LPS2IECDirections,DirectionType); PointType GetCOV(); + PointType GetOriginWOffset(); protected: @@ -210,7 +235,8 @@ protected: m_Size; PointType - m_OriginLPS; + m_OriginLPS, + m_ImportOffset; DirectionType m_ImageDirections, @@ -239,6 +265,8 @@ public: typedef itk::Object Superclass; typedef itk::SmartPointer Pointer; typedef itk::Point PointType; + typedef itk::Vector SpacingType; + typedef itk::Size<3> SizeType; /** Method for creation through the object factory. */ @@ -262,11 +290,23 @@ public: itkSetMacro(IntensityThreshold, double); itkGetMacro(IntensityThreshold, double); - itkSetMacro(ProjectionOriginLPS,PointType); - itkGetMacro(ProjectionOriginLPS,PointType); + itkSetMacro(DRT1Size, SizeType); + itkGetMacro(DRT1Size, SizeType); - itkSetMacro(ProjectionOriginLPSZero,PointType); - itkGetMacro(ProjectionOriginLPSZero,PointType); + itkSetMacro(DRT2Size, SizeType); + itkGetMacro(DRT2Size, SizeType); + + itkSetMacro(DRT1Spacing, SpacingType); + itkGetMacro(DRT1Spacing, SpacingType); + + itkSetMacro(DRT2Spacing, SpacingType); + itkGetMacro(DRT2Spacing, SpacingType); + + itkSetMacro(IECS2IECScannerT, PointType); + itkGetMacro(IECS2IECScannerT, PointType); + + itkSetMacro(IECS2IECScannerR, PointType); + itkGetMacro(IECS2IECScannerR, PointType); itkSetMacro(ProjectionCenter, PointType); itkGetMacro(ProjectionCenter, PointType); @@ -279,9 +319,20 @@ protected: m_SCD, m_IntensityThreshold; + + SizeType + m_DRT1Size, + m_DRT2Size; + + SpacingType + m_DRT1Spacing, + m_DRT2Spacing; + + PointType + m_IECS2IECScannerT, + m_IECS2IECScannerR; + PointType - m_ProjectionOriginLPS, - m_ProjectionOriginLPSZero, /* center of projection in an IEC reference at * Patient Origin of fixed images. Positioning scanner */ m_ProjectionCenter; @@ -301,6 +352,55 @@ private: }; + + +class RTGeometryMetaInformation : + public itk::Object{ + +public: + /** standard typedefs **/ + typedef RTGeometryMetaInformation Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::Point PointType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(RTGeometryMetaInformation, itk::Object); + + /** object information streaming **/ + void PrintSelf(std::ostream& os, itk::Indent indent) const; + + + itkSetMacro(IsocenterLPS,PointType); + itkGetMacro(IsocenterLPS,PointType); + + itkSetMacro(IsocenterIECS,PointType); + itkGetMacro(IsocenterIECS,PointType); + + + +protected: + + PointType + m_IsocenterLPS, + m_IsocenterIECS; + + /** Default Constructor **/ + RTGeometryMetaInformation (); + /** Default Destructor **/ + virtual ~RTGeometryMetaInformation (); + +private: + /** purposely not implemented **/ + RTGeometryMetaInformation (const Self&); + + /** purposely not implemented **/ + void operator=(const Self&); +}; + } #endif diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 50e98af..15e3ec0 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -98,9 +98,9 @@ itkImageProcessor::itkImageProcessor() interpolator2 = InterpolatorType::New(); - image1res[0]=image2res[1]= 2.; - image1Size[0] = image1Size[1] = 512; - image2Size[0] = image2Size[1] = 512; +// image1res[0]=image2res[1]= 2.; +// image1Size[0] = image1Size[1] = 512; +// image2Size[0] = image2Size[1] = 512; TZero[0]=TZero[1]=TZero[2]=0.; @@ -137,10 +137,10 @@ itkImageProcessor::itkImageProcessor() m_3DInputChangeInformationToZero = ChangeInformationFilterType::New(); - IEC2DCMMapT.Fill(0.); - IEC2DCMMapR.Fill(0.); - rtIsocenterLPS.Fill(0.); - rtCouchOffset.Fill(0.); +// IEC2DCMMapT.Fill(0.); +// IEC2DCMMapR.Fill(0.); +// rtIsocenterLPS.Fill(0.); +// rtCouchOffset.Fill(0.); //ProjectionCenterFixedAxes.Fill(0.); @@ -151,18 +151,50 @@ itkImageProcessor::itkImageProcessor() + /* Set to NULL the metainfo that are filled in on load */ m_CTMetaInfo = NULL; + m_TImage1MetaInfo = NULL; m_TImage2MetaInfo = NULL; + m_DRTImage1MetaInfo = NULL; + m_DRTImage2MetaInfo = NULL; + + m_RTMetaInfo = NULL; + + /* Initialiser the projection geoemtry with defaults */ m_DRTGeometryMetaInfo = DRTProjectionGeometryImageMetaInformation::New(); m_DRTGeometryMetaInfo->SetSCD(570.); m_DRTGeometryMetaInfo->SetProjectionAngle1IEC(180.); m_DRTGeometryMetaInfo->SetProjectionAngle2IEC(90.); + m_DRTGeometryMetaInfo->SetIntensityThreshold(-300.); + ImageType3D::PointType + Point3D; + Point3D[0]=0.; + Point3D[1]=0.; + Point3D[2]=175.; + m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D); + ImageType3D::SizeType + ImageSize; + ImageSize[0]=512; + ImageSize[1]=512; + ImageSize[2]=1; + m_DRTGeometryMetaInfo->SetDRT1Size(ImageSize); + m_DRTGeometryMetaInfo->SetDRT2Size(ImageSize); + ImageType3D::SpacingType + ImageRes; + ImageRes[0] = 1.; + ImageRes[1] = 1.; + ImageRes[2] = 1.; + m_DRTGeometryMetaInfo->SetDRT1Spacing(ImageRes); + m_DRTGeometryMetaInfo->SetDRT2Spacing(ImageRes); + Point3D[0]=0.; + Point3D[1]=0.; + Point3D[2]=0.; + m_DRTGeometryMetaInfo->SetIECS2IECScannerT(Point3D); + m_DRTGeometryMetaInfo->SetIECS2IECScannerR(Point3D); - m_DRTImage1MetaInfo = DRTImageMetaInformation::New(); - m_DRTImage2MetaInfo = DRTImageMetaInformation::New(); @@ -201,34 +233,49 @@ void itkImageProcessor::SetCustom_ImportTransform(double dTx, double dRy, double dRz) { - IEC2DCMMapT[0] = dTx; - IEC2DCMMapT[1] = dTy; - IEC2DCMMapT[2] = dTz; - IEC2DCMMapR[0] = dRx; - IEC2DCMMapR[1] = dRy; - IEC2DCMMapR[2] = dRz; +// IEC2DCMMapT[0] = dTx; +// IEC2DCMMapT[1] = dTy; +// IEC2DCMMapT[2] = dTz; +// IEC2DCMMapR[0] = dRx; +// IEC2DCMMapR[1] = dRy; +// IEC2DCMMapR[2] = dRz; customized_ImportTransform = true; + + ImageType3D::PointType + Punto; + Punto[0] = dTx; + Punto[1] = dTy; + Punto[2] = dTz; + + m_DRTGeometryMetaInfo->SetIECS2IECScannerT(Punto); + + Punto[0] = dRx; + Punto[1] = dRy; + Punto[2] = dRz; + + m_DRTGeometryMetaInfo->SetIECS2IECScannerR(Punto); + } -void itkImageProcessor::SetCustom_RTCouchSetup_IEC(double dLAT, - double dVRT, - double dLNG) -{ - rtCouchOffset[0] = dLAT; - rtCouchOffset[1] = dVRT; - rtCouchOffset[2] = dLNG; - customized_RTCouchSetup = true; -} +//void itkImageProcessor::SetCustom_RTCouchSetup_IEC(double dLAT, +// double dVRT, +// double dLNG) +//{ +// rtCouchOffset[0] = dLAT; +// rtCouchOffset[1] = dVRT; +// rtCouchOffset[2] = dLNG; +// customized_RTCouchSetup = true; +//} -void itkImageProcessor::SetCustom_RTIsocenter_LPS(double dX, - double dY, - double dZ) -{ - rtIsocenterLPS[0] = dX; - rtIsocenterLPS[1] = dY; - rtIsocenterLPS[2] = dZ; - customized_RTIsocenter = true; -} +//void itkImageProcessor::SetCustom_RTIsocenter_LPS(double dX, +// double dY, +// double dZ) +//{ +// rtIsocenterLPS[0] = dX; +// rtIsocenterLPS[1] = dY; +// rtIsocenterLPS[2] = dZ; +// customized_RTIsocenter = true; +//} void itkImageProcessor::SetCustom_ProjectionCenterFixedAxes_IEC(double dX, double dY, @@ -249,20 +296,53 @@ void itkImageProcessor::SetCustom_ProjectionCenterFixedAxes_IEC(double dX, void itkImageProcessor::SetCustom_2Dres(double nX1,double nY1,double nX2,double nY2) { - image1res[0] = nX1; - image1res[1] = nY1; - image2res[0] = nX2; - image2res[1] = nY2; +// image1res[0] = nX1; +// image1res[1] = nY1; +// image2res[0] = nX2; +// image2res[1] = nY2; customized_2DRES = true; + + if(m_DRTGeometryMetaInfo == NULL) { + // todo + } + + ImageType3D::SpacingType Spacing; + Spacing [0] = nX1; + Spacing [1] = nY1; + Spacing [1] = 1.; + m_DRTGeometryMetaInfo->SetDRT1Spacing(Spacing); + + Spacing [0] = nX2; + Spacing [1] = nY2; + Spacing [1] = 1.; + m_DRTGeometryMetaInfo->SetDRT2Spacing(Spacing); + //TODO UPDATE TO FOLLOW } void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2) { - image1Size[0] = nX1; - image1Size[1] = nY1; - image2Size[0] = nX2; - image2Size[1] = nY2; +// image1Size[0] = nX1; +// image1Size[1] = nY1; +// image2Size[0] = nX2; +// image2Size[1] = nY2; customized_2DSIZE = true; + + if(m_DRTGeometryMetaInfo == NULL) { + // todo + } + + ImageType3D::SizeType Size; + Size [0] = nX1; + Size [1] = nY1; + Size [1] = 1.; + m_DRTGeometryMetaInfo->SetDRT1Size(Size); + + Size [0] = nX2; + Size [1] = nY2; + Size [1] = 1.; + m_DRTGeometryMetaInfo->SetDRT2Size(Size); + + //TODO UPDATE TO FOLLOW } @@ -393,33 +473,61 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) InternalImageType::Pointer m_InputImage = m_VolumeSourceDupli->GetOutput(); + if(m_CTMetaInfo != NULL){ // TODO UNLOAD } + if(m_DRTImage1MetaInfo != NULL){ + // TODO UNLOAD + } + + if(m_DRTImage2MetaInfo != NULL){ + // TODO UNLOAD + } + + if(m_RTMetaInfo != NULL){ + // TODO UNLOAD + } + + + /* copy useful meta information into the CT container */ m_CTMetaInfo = CTVolumeImageMetaInformation::New(); + m_CTMetaInfo->SetPatientOrientation(m_PatOrientation); + m_CTMetaInfo->SetSpacing(m_InputImage->GetSpacing()); m_CTMetaInfo->SetSize( m_InputImage->GetBufferedRegion().GetSize() ); m_CTMetaInfo->SetOriginLPS(m_InputImage->GetOrigin()); - m_CTMetaInfo->SetSpacing(m_InputImage->GetSpacing()); m_CTMetaInfo->SetImageDirections(m_InputImage->GetDirection()); - m_CTMetaInfo->SetPatientOrientation(m_PatOrientation); + ImageType3D::PointType + PointOffset; + PointOffset.Fill(0.); + m_CTMetaInfo->SetImportOffset(PointOffset); + /* initialise DRT meta */ + m_DRTImage1MetaInfo = DRTImageMetaInformation::New(); + m_DRTImage1MetaInfo->SetProjectionAngleLPS( + CalcProjectionAngleLPS( + m_CTMetaInfo->GetPatientOrientation(), + m_DRTGeometryMetaInfo->GetProjectionAngle1IEC()) + ); + m_DRTImage1MetaInfo->SetSCD( + m_DRTGeometryMetaInfo->GetSCD()); /* Calculate projection angle IEC to LPS */ - m_DRTImage1MetaInfo->SetProjectionAngleLPS( - CalcProjectionAngleLPS( - m_CTMetaInfo->GetPatientOrientation(), - m_DRTGeometryMetaInfo->GetProjectionAngle1IEC()) - ); - + m_DRTImage2MetaInfo = DRTImageMetaInformation::New(); m_DRTImage2MetaInfo->SetProjectionAngleLPS( CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) - ); + ); + m_DRTImage2MetaInfo->SetSCD( + m_DRTGeometryMetaInfo->GetSCD()); + + + this->UpdateProjectionGeometryMeta(); @@ -443,6 +551,9 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) m_3DInputChangeInformationToZero->GetOutput(); + + + /* TODO: Here calculate COV * m_3DProjectionOriginLPSZero * m_3DProjectionOriginLPS */ @@ -453,116 +564,185 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) } +itkImageProcessor::ImageType3D::PointType +itkImageProcessor::CalculateDRTOrigin( + ImageType3D::SizeType m_ImageSize, + ImageType3D::SpacingType m_ImageResolution, + double dSCD + ){ -void itkImageProcessor::ApplyVolumeImportTransform(){ + ImageType3D::PointType Origin; + double o2Dx, o2Dy; + o2Dx = ((double)m_ImageSize [0] - 1.) / 2.; + o2Dy = ((double)m_ImageSize [1] - 1.) / 2.; + // Compute the origin (in mm) of the 2D Image + Origin[0] = -m_ImageResolution [0] * o2Dx; + Origin[1] = -m_ImageResolution [1] * o2Dy; + Origin[2] = -dSCD ; + return + Origin; +} +itkImageProcessor::ImageType3D::PointType +itkImageProcessor::CalculateProjectionCenterLPS( + ImageType3D::PointType m_ImportOffset, + ImageType3D::DirectionType m_VolumeLPS2IEC, + ImageType3D::PointType m_VolumeCOVLPS, + ImageType3D::PointType m_VolumeOriginLPS, + ImageType3D::PointType m_ProjectionCenter, + bool bZero) +{ - //ImageType3D::PointType ImportOffset; - - if(customized_ImportTransform && - customized_RTCouchSetup && - customized_RTIsocenter ){ - - ImportOffset= - this->CalcImportVolumeOffset(rtCouchOffset, - m_CTMetaInfo->GetLPS2IECDirections(), - rtIsocenterLPS, - IEC2DCMMapT); - std::cout<<"ImportOffset> "<GetOutput()->GetOrigin(); - Origin = Origin - ImportOffset ; + Origin = m_VolumeOriginLPS; + Origin = Origin - m_ImportOffset; - ChangeInformationFilterType::Pointer - m_3DInputChangeInformation = - ChangeInformationFilterType::New(); - m_3DInputChangeInformation->SetInput( - m_VolumeSourceDupli->GetOutput()); - m_3DInputChangeInformation->SetOutputOrigin(Origin ) ;//NewOrigin); - m_3DInputChangeInformation->ChangeOriginOn(); - m_3DInputChangeInformation->UpdateOutputInformation(); - m_3DInputChangeInformation->Update(); - - /* End of origin change */ - - InternalImageType::Pointer m_Image = - m_3DInputChangeInformation->GetOutput(); - - // TODO: move this view origin to new meta variable. - m_CTMetaInfo->SetOriginLPS(m_Image->GetOrigin()); - - /* Calculate center of projection in patient coordinates */ InternalImageType::DirectionType IECtoLPS_Directions; - IECtoLPS_Directions = - m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); - // LPStoIEC_Directions.GetTranspose(); + IECtoLPS_Directions = m_VolumeLPS2IEC.GetTranspose(); - if(customized_ProjectionCenter) + ImageType3D::PointType m_ProjectionOriginLPS; + m_ProjectionOriginLPS = + IECtoLPS_Directions * m_ProjectionCenter; + /* in longitudinal direction (Sup-Inf), we put it in the + * middle of the volume */ + m_ProjectionOriginLPS[2] = m_VolumeCOVLPS[2] - m_ImportOffset[2]; + + if(bZero) { - m_3DProjectionOriginLPS = - IECtoLPS_Directions * m_DRTGeometryMetaInfo->GetProjectionCenter(); - /* in longitudinal direction (Sup-Inf), we put it in the - * middle of the volume */ - m_3DProjectionOriginLPS[2] = m_CTMetaInfo->GetCOV()[2]; - + return + m_ProjectionOriginLPS - Origin; } else { - m_3DProjectionOriginLPS = m_CTMetaInfo->GetCOV(); + return + m_ProjectionOriginLPS; } - - - /* calculate also the same center with zero origin - * This is required to bring back the DRT into LPS */ - m_3DProjectionOriginLPSZero = - m_3DProjectionOriginLPS - m_CTMetaInfo->GetOriginLPS() ; - - /* We should know everything we need to calculate the - * origin and direction of projection images */ - - //std::cout<< "m_3DProjectionOriginLPS "<GetCOV() <CalcImportVolumeOffset(rtCouchOffset, +// m_CTMetaInfo->GetLPS2IECDirections(), +// rtIsocenterLPS, +// IEC2DCMMapT); +// std::cout<<"ImportOffset> "<GetOutput()->GetOrigin(); +// Origin = Origin - ImportOffset ; + +// ChangeInformationFilterType::Pointer +// m_3DInputChangeInformation = +// ChangeInformationFilterType::New(); +// m_3DInputChangeInformation->SetInput( +// m_VolumeSourceDupli->GetOutput()); +// m_3DInputChangeInformation->SetOutputOrigin(Origin ) ;//NewOrigin); +// m_3DInputChangeInformation->ChangeOriginOn(); +// m_3DInputChangeInformation->UpdateOutputInformation(); +// m_3DInputChangeInformation->Update(); + +// /* End of origin change */ + +// InternalImageType::Pointer m_Image = +// m_3DInputChangeInformation->GetOutput(); + +// // TODO: move this view origin to new meta variable. +// m_CTMetaInfo->SetOriginLPS(m_Image->GetOrigin()); + +// /* Calculate center of projection in patient coordinates */ + +// InternalImageType::DirectionType IECtoLPS_Directions; +// IECtoLPS_Directions = +// m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); +// // LPStoIEC_Directions.GetTranspose(); + +// if(customized_ProjectionCenter) +// { +// m_3DProjectionOriginLPS = +// IECtoLPS_Directions * m_DRTGeometryMetaInfo->GetProjectionCenter(); +// /* in longitudinal direction (Sup-Inf), we put it in the +// * middle of the volume */ +// m_3DProjectionOriginLPS[2] = m_CTMetaInfo->GetCOV()[2]; + +// } else { +// m_3DProjectionOriginLPS = m_CTMetaInfo->GetCOV(); +// } + + +// /* calculate also the same center with zero origin +// * This is required to bring back the DRT into LPS */ +// m_3DProjectionOriginLPSZero = +// m_3DProjectionOriginLPS - m_CTMetaInfo->GetOriginLPS() ; + +// /* We should know everything we need to calculate the +// * origin and direction of projection images */ + +// //std::cout<< "m_3DProjectionOriginLPS "<GetCOV() <GetProjectionOriginLPS()<GetProjectionOriginLPSZero()<GetProjectionOriginLPS()<GetProjectionOriginLPSZero()< itkImageProcessor::GetLPStoProjectionGeoLPSOffset() { - ConversionOffset[0] = - - m_3DProjectionOriginLPS[0] + m_3DProjectionOriginLPSZero[0]; - ConversionOffset[1] = - - m_3DProjectionOriginLPS[1] + m_3DProjectionOriginLPSZero[1]; - ConversionOffset[2] = - - m_3DProjectionOriginLPS[2] + m_3DProjectionOriginLPSZero[2]; + + std::vector vOffset; + + vOffset.push_back( + m_DRTImage1MetaInfo->GetLPStoProjectionGeoLPSOffset()[0] + ); + vOffset.push_back( + m_DRTImage1MetaInfo->GetLPStoProjectionGeoLPSOffset()[1] + ); + vOffset.push_back( + m_DRTImage1MetaInfo->GetLPStoProjectionGeoLPSOffset()[2] + ); return - ConversionOffset.GetDataPointer(); + vOffset; + } double @@ -853,9 +1033,8 @@ int itkImageProcessor::load2D(const char * pcFName){ currImgOrient); m_TImageMeta->SetProjectionOrientation( currProjOrientation); - // memcpy(m_ImageIntensity,dIntensityWindow, 2*sizeof(double)); - m_Duplicator->SetInputImage(caster2DInput->GetOutput() );//filter->GetOutput()); + m_Duplicator->SetInputImage(caster2DInput->GetOutput() ); m_Duplicator->Update(); @@ -946,14 +1125,20 @@ void itkImageProcessor::InitializeProjector() dtr * RZero[1], dtr * RZero[2]); - transform->SetCenter(m_3DProjectionOriginLPSZero); + if(m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() != + m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ) { + // TODO + } + + transform->SetCenter( + m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() ); // 2D Image 1 interpolator1->SetProjectionAngle( dtr * m_DRTImage1MetaInfo->GetProjectionAngleLPS() ); interpolator1->SetFocalPointToIsocenterDistance( - m_DRTGeometryMetaInfo->GetSCD()); + m_DRTImage1MetaInfo->GetSCD()); interpolator1->SetThreshold( m_DRTGeometryMetaInfo->GetIntensityThreshold() ); @@ -967,7 +1152,7 @@ void itkImageProcessor::InitializeProjector() dtr * m_DRTImage2MetaInfo->GetProjectionAngleLPS() ); interpolator2->SetFocalPointToIsocenterDistance( - m_DRTGeometryMetaInfo->GetSCD()); + m_DRTImage2MetaInfo->GetSCD()); interpolator2->SetThreshold( m_DRTGeometryMetaInfo->GetIntensityThreshold() ); @@ -1001,6 +1186,153 @@ itkImageProcessor::CalcDRTImageDirections( } +void itkImageProcessor::loadRTPlanInfo( + double dIsoX, double dIsoY, double dIsoZ, + double dLAT, double dVRT ,double dLNG){ + + if(m_RTMetaInfo != NULL){ + //TODO + } + + m_RTMetaInfo = RTGeometryMetaInformation::New(); + + ImageType3D::PointType + Point; + Point[0] = dIsoX; + Point[1] = dIsoY; + Point[2] = dIsoZ; + + m_RTMetaInfo->SetIsocenterLPS(Point); + + Point[0] = dLAT; + Point[1] = dVRT; + Point[2] = dLNG; + + m_RTMetaInfo->SetIsocenterIECS(Point); + + m_CTMetaInfo->SetImportOffset( + CalcImportVolumeOffset( + m_RTMetaInfo->GetIsocenterIECS(), + m_CTMetaInfo->GetLPS2IECDirections(), + m_RTMetaInfo->GetIsocenterLPS(), + m_DRTGeometryMetaInfo->GetIECS2IECScannerT() + ) + ); + + std::cout<< "m_CTMetaInfo->GetImportOffset() " + <GetImportOffset() <UpdateProjectionGeometryMeta(); + + +} + +void itkImageProcessor::UpdateProjectionGeometryMeta(){ + + if(m_CTMetaInfo == NULL){ + //TODO. + } + + m_DRTImage1MetaInfo->SetProjectionOriginLPS( + CalculateProjectionCenterLPS( + m_CTMetaInfo->GetImportOffset(), + m_CTMetaInfo->GetLPS2IECDirections(), + m_CTMetaInfo->GetCOV(), + m_CTMetaInfo->GetOriginLPS(), + m_DRTGeometryMetaInfo->GetProjectionCenter(), + false + ) + ); + m_DRTImage1MetaInfo->SetProjectionOriginLPSZero( + CalculateProjectionCenterLPS( + m_CTMetaInfo->GetImportOffset(), + m_CTMetaInfo->GetLPS2IECDirections(), + m_CTMetaInfo->GetCOV(), + m_CTMetaInfo->GetOriginLPS(), + m_DRTGeometryMetaInfo->GetProjectionCenter(), + true + ) + ); + //TODO: Calculate overlapping region + m_DRTImage1MetaInfo->SetSize(m_DRTGeometryMetaInfo->GetDRT1Size()); + m_DRTImage1MetaInfo->SetSpacing(m_DRTGeometryMetaInfo->GetDRT1Spacing()); + m_DRTImage1MetaInfo->SetOrigin( + CalculateDRTOrigin( + m_DRTImage1MetaInfo->GetSize(), + m_DRTImage1MetaInfo->GetSpacing(), + m_DRTImage1MetaInfo->GetSCD() + ) + ); + m_DRTImage1MetaInfo->SetOriginLPS( + CalcDRTImageOrigin( + m_DRTImage1MetaInfo->GetOrigin(), + m_DRTImage1MetaInfo->GetCOV(), + m_DRTImage1MetaInfo->GetProjectionOriginLPS(), + m_DRTImage1MetaInfo->GetProjectionAngleLPS(), + Std_DRT2LPS + ) + ); + m_DRTImage1MetaInfo->SetImageDirectionsLPS( + CalcDRTImageDirections( + m_DRTImage1MetaInfo->GetProjectionAngleLPS(), + Std_DRT2LPS) + ); + + + m_DRTImage2MetaInfo->SetProjectionOriginLPS( + CalculateProjectionCenterLPS( + m_CTMetaInfo->GetImportOffset(), + m_CTMetaInfo->GetLPS2IECDirections(), + m_CTMetaInfo->GetCOV(), + m_CTMetaInfo->GetOriginLPS(), + m_DRTGeometryMetaInfo->GetProjectionCenter(), + false + ) + ); + m_DRTImage2MetaInfo->SetProjectionOriginLPSZero( + CalculateProjectionCenterLPS( + m_CTMetaInfo->GetImportOffset(), + m_CTMetaInfo->GetLPS2IECDirections(), + m_CTMetaInfo->GetCOV(), + m_CTMetaInfo->GetOriginLPS(), + m_DRTGeometryMetaInfo->GetProjectionCenter(), + true + ) + ); + //TODO: Calculate overlapping region + m_DRTImage2MetaInfo->SetSize(m_DRTGeometryMetaInfo->GetDRT2Size()); + m_DRTImage2MetaInfo->SetSpacing(m_DRTGeometryMetaInfo->GetDRT2Spacing()); + m_DRTImage2MetaInfo->SetOrigin( + CalculateDRTOrigin( + m_DRTImage2MetaInfo->GetSize(), + m_DRTImage2MetaInfo->GetSpacing(), + m_DRTImage2MetaInfo->GetSCD() + ) + ); + m_DRTImage2MetaInfo->SetOriginLPS( + CalcDRTImageOrigin( + m_DRTImage2MetaInfo->GetOrigin(), + m_DRTImage2MetaInfo->GetCOV(), + m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + m_DRTImage2MetaInfo->GetProjectionAngleLPS(), + Std_DRT2LPS + ) + ); + m_DRTImage2MetaInfo->SetImageDirectionsLPS( + CalcDRTImageDirections( + m_DRTImage2MetaInfo->GetProjectionAngleLPS(), + Std_DRT2LPS) + ); + + std::cout<< "Projection Origin LPS" <GetProjectionOriginLPS() <GetOriginLPS() <GetImageDirectionsLPS() <GetProjectionOriginLPS() <GetOriginLPS() <GetImageDirectionsLPS() <CalcDRTImageDirections(dAngle, stdDRT2LPS); + + ImageType3D::PointType NewOrigin = + m_ProjectionCenterLPS + DRT2LPS * (m_DRTOrigin - m_DRTCOV); + + return + NewOrigin; +} + void itkImageProcessor::GetProjectionImages(){ @@ -1070,8 +1421,15 @@ void itkImageProcessor::GetProjectionImages(){ finalTransform->SetRotation(dtr * RZero[0], dtr * RZero[1], dtr * RZero[2]); // The transform is determined by the parameters and the rotation center. + + if(m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() != + m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ) { + // TODO + } + finalTransform->SetCenter( - m_3DProjectionOriginLPSZero); + m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()); + using ResampleFilterType = itk::ResampleImageFilter; @@ -1090,31 +1448,9 @@ void itkImageProcessor::GetProjectionImages(){ interpolator1->Initialize(); resampleFilter1->SetInterpolator(interpolator1); - // The output 2D projection image has the same image size, origin, and the pixel spacing as - // those of the input 2D image. - InternalImageType::SizeType size; - double spacing[Dimension]; - double o2Dx, o2Dy, origin[Dimension]; - - - size[0] = image1Size[0]; // number of pixels along X of the 2D DRR image - size[1] = image1Size[1]; // number of pixels along X of the 2D DRR image - size[2] = 1; // only one slice - spacing[0] = image1res[0]; // pixel spacing along X of the 2D DRR image [mm] - spacing[1] = image1res[1]; // pixel spacing along Y of the 2D DRR image [mm] - spacing[2] = 1.0; // slice thickness of the 2D DRR image [mm] - - o2Dx = ((double)image1Size[0] - 1.) / 2.; - o2Dy = ((double)image1Size[1] - 1.) / 2.; - // Compute the origin (in mm) of the 2D Image - origin[0] = -image1res[0] * o2Dx; - origin[1] = -image1res[1] * o2Dy; - origin[2] = -m_DRTGeometryMetaInfo->GetSCD(); - - resampleFilter1->SetSize(size); - resampleFilter1->SetOutputSpacing(spacing); - resampleFilter1->SetOutputOrigin(origin); - + resampleFilter1->SetSize(m_DRTImage1MetaInfo->GetSize()); + resampleFilter1->SetOutputSpacing(m_DRTImage1MetaInfo->GetSpacing()); + resampleFilter1->SetOutputOrigin(m_DRTImage1MetaInfo->GetOrigin()); // Do the same thing for the output image 2. ResampleFilterType::Pointer resampleFilter2 = ResampleFilterType::New(); @@ -1131,26 +1467,9 @@ void itkImageProcessor::GetProjectionImages(){ interpolator2->Initialize(); resampleFilter2->SetInterpolator(interpolator2); - - - size[0] = image2Size[0]; // number of pixels along X of the 2D DRR image - size[1] = image2Size[1]; // number of pixels along X of the 2D DRR image - size[2] = 1; // only one slice - spacing[0] = image2res[0]; // pixel spacing along X of the 2D DRR image [mm] - spacing[1] = image2res[1]; // pixel spacing along Y of the 2D DRR image [mm] - spacing[2] = 1.0; // slice thickness of the 2D DRR image [mm] - - o2Dx = ((double)image2Size[0] - 1.) / 2.; - o2Dy = ((double)image2Size[1] - 1.) / 2.; - // Compute the origin (in mm) of the 2D Image - origin[0] = -image2res[0] * o2Dx; - origin[1] = -image2res[1] * o2Dy; - origin[2] = -m_DRTGeometryMetaInfo->GetSCD(); - - resampleFilter2->SetSize(size); - resampleFilter2->SetOutputSpacing(spacing); - resampleFilter2->SetOutputOrigin(origin); - + resampleFilter2->SetSize(m_DRTImage2MetaInfo->GetSize()); + resampleFilter2->SetOutputSpacing(m_DRTImage2MetaInfo->GetSpacing()); + resampleFilter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOrigin()); resampleFilter1->Update(); @@ -1173,26 +1492,11 @@ void itkImageProcessor::GetProjectionImages(){ /* First of all we need the center of the Projection images in its reference frame */ - - ImageType3D::PointType image3DOrigin = m_3DProjectionOriginLPS; - - resampleFilter1->Update(); - - ImageType3D::PointType NewOrigin = - this->CalcDRTImageOrigin(resampleFilter1->GetOutput(), - m_3DProjectionOriginLPS, - m_DRTImage1MetaInfo->GetProjectionAngleLPS(), - Std_DRT2LPS); - InternalImageType::DirectionType DRT2LPS1 = - this->CalcDRTImageDirections( - m_DRTImage1MetaInfo->GetProjectionAngleLPS(), - Std_DRT2LPS); - filter1->SetInput( resampleFilter1->GetOutput() ); - filter1->SetOutputDirection(DRT2LPS1 );//IECtoLPS_Directions); + filter1->SetOutputDirection(m_DRTImage1MetaInfo->GetImageDirectionsLPS() );//IECtoLPS_Directions); filter1->ChangeDirectionOn(); - filter1->SetOutputOrigin( NewOrigin ); + filter1->SetOutputOrigin( m_DRTImage1MetaInfo->GetOriginLPS() ); filter1->ChangeOriginOn(); filter1->UpdateOutputInformation(); @@ -1200,24 +1504,10 @@ void itkImageProcessor::GetProjectionImages(){ /* Again.. */ resampleFilter2->Update(); - - NewOrigin = - this->CalcDRTImageOrigin(resampleFilter2->GetOutput(), - m_3DProjectionOriginLPS, - m_DRTImage2MetaInfo->GetProjectionAngleLPS(), - Std_DRT2LPS); - - DRT2LPS1 = - this->CalcDRTImageDirections( - m_DRTImage2MetaInfo->GetProjectionAngleLPS(), - Std_DRT2LPS); - - - filter2->SetInput( resampleFilter2 ->GetOutput() ); - filter2->SetOutputDirection(DRT2LPS1 ); + filter2->SetOutputDirection(m_DRTImage2MetaInfo->GetImageDirectionsLPS() ); filter2->ChangeDirectionOn(); - filter2->SetOutputOrigin(NewOrigin ); + filter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOriginLPS() ); filter2->ChangeOriginOn(); filter2->UpdateOutputInformation(); diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index 2be363e..c660bfc 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -75,6 +75,7 @@ public: int load3DSerie(const char* ); int load2D(const char *); + void loadRTPlanInfo(double, double, double, double, double ,double); /** Projection angles - Gantry angle IEC */ void SetProjectionAngle1IEC(double); @@ -96,8 +97,8 @@ public: /** Custom settings of the projection images */ void SetCustom_2Dres(double,double,double,double); void SetCustom_2Dsize(int,int,int,int); - void SetCustom_RTIsocenter_LPS(double, double, double); - void SetCustom_RTCouchSetup_IEC(double, double, double); + //void SetCustom_RTIsocenter_LPS(double, double, double); + //void SetCustom_RTCouchSetup_IEC(double, double, double); void SetCustom_ImportTransform(double, double, double, double, double, double); void SetCustom_ProjectionCenterFixedAxes_IEC(double,double,double); @@ -110,7 +111,7 @@ public: void InitializeProjector(); /** Get Projection origin in LPS coordinates*/ - const double* GetLPStoProjectionGeoLPSOffset(); + const std::vector GetLPStoProjectionGeoLPSOffset(); /** Get Projection origin in LPS coordinates*/ const double* GetRTImportOffset(); @@ -119,7 +120,10 @@ public: /** Apply transform to CT volume. * This moves the volume based on RT Plan info * if those are available. */ - void ApplyVolumeImportTransform(); + //void ApplyVolumeImportTransform(); + + + /** Get the Localizer intensity window and * center for rendering */ @@ -191,42 +195,61 @@ private: using ITKtoVTKFilterType = itk::ImageToVTKImageFilter; -ImageType3D::PointType ImportOffset; + ImageType3D::PointType + CalculateDRTOrigin( + ImageType3D::SizeType m_ImageSize, + ImageType3D::SpacingType m_ImageResolution, + double dSCD + ); + + ImageType3D::PointType + CalculateProjectionCenterLPS( + ImageType3D::PointType m_ImportOffset, + ImageType3D::DirectionType m_VolumeLPS2IEC, + ImageType3D::PointType m_VolumeCOVLPS, + ImageType3D::PointType m_VolumeOriginLPS, + ImageType3D::PointType m_ProjectionCenter, + bool bZero); + + + ImageType3D::PointType ImportOffset; ImageType3D :: PointType - ConversionOffset; + ConversionOffset; TransformType::Pointer - transform; + transform; InterpolatorType::Pointer - interpolator1, - interpolator2; + interpolator1, + interpolator2; ITKtoVTKFilterType::Pointer - toVTK2D1, - toVTK2D2, - toVTKLocalizer1, - toVTKLocalizer2; + toVTK2D1, + toVTK2D2, + toVTKLocalizer1, + toVTKLocalizer2; InternalImageType::Pointer - image3DIn, - imageDRT1In, - imageDRT2In; + image3DIn, + imageDRT1In, + imageDRT2In; DuplicatorType::Pointer - m_LATSourceDupli, - m_PASourceDupli, - m_VolumeSourceDupli; + m_LATSourceDupli, + m_PASourceDupli, + m_VolumeSourceDupli; ChangeInformationFilterType::Pointer - m_3DInputChangeInformationToZero; + m_3DInputChangeInformationToZero; - ImageType3D :: PointType - m_3DProjectionOriginLPS, - m_3DProjectionOriginLPSZero; +// ImageType3D :: PointType +// m_3DProjectionOriginLPS, +// m_3DProjectionOriginLPSZero; + void UpdateProjectionGeometryMeta(); + /* The following functions do not rely on class variables, * Only input variables are used... */ @@ -241,6 +264,16 @@ ImageType3D::PointType ImportOffset; InternalImageType::DirectionType stdDRT2LPS ); + ImageType3D::PointType + CalcDRTImageOrigin( + ImageType3D::PointType m_DRTOrigin, + ImageType3D::PointType m_DRTCOV, + ImageType3D::PointType m_ProjectionCenterLPS, + double dAngle, + InternalImageType::DirectionType stdDRT2LPS + ); + + double CalcProjectionAngleLPS( tPatOrientation pOrient, @@ -253,8 +286,8 @@ ImageType3D::PointType ImportOffset; ImageType3D::PointType rtIsocenterLPS, ImageType3D::PointType IEC2DCMMapT); - double image1res[2], image2res[2]; - int image1Size[2], image2Size[2]; +// double image1res[2], image2res[2]; +// int image1Size[2], image2Size[2]; double TZero[3], RZero[3]; @@ -276,23 +309,23 @@ ImageType3D::PointType ImportOffset; /* Transform between IEC Support and * IEC Scanner frame of reference */ - ImageType3D::PointType - IEC2DCMMapT, - IEC2DCMMapR; +// ImageType3D::PointType +// IEC2DCMMapT, +// IEC2DCMMapR; /* RT Plan isocenter in LPS coordinates */ - ImageType3D::PointType - rtIsocenterLPS; +// ImageType3D::PointType +// rtIsocenterLPS; /* RT Plan couch SETUP offset corresponding to * the RT Plan isocenter */ - ImageType3D::PointType - rtCouchOffset; +// ImageType3D::PointType +// rtCouchOffset; InternalImageType::DirectionType - Std_DRT2LPS; + Std_DRT2LPS; @@ -301,6 +334,10 @@ ImageType3D::PointType ImportOffset; * m_Projection2VTKTransform; + + /** + * Many meta containers + */ CTVolumeImageMetaInformation::Pointer m_CTMetaInfo; @@ -315,6 +352,9 @@ ImageType3D::PointType ImportOffset; m_TImage1MetaInfo, m_TImage2MetaInfo; + RTGeometryMetaInformation::Pointer + m_RTMetaInfo; + };