diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp index 6f8a6d2..8368cc7 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp @@ -80,7 +80,7 @@ DRTImageMetaInformation(){ this->m_Spacing.Fill(0.); - this->m_Origin.Fill(0.); + //this->m_Origin.Fill(0.); this->m_ProjectionOriginLPS.Fill(0.); @@ -111,6 +111,98 @@ DRTImageMetaInformation } +DRTImageMetaInformation::PointType +DRTImageMetaInformation::GetOrigin() +{ + + //PointType Origin; + PointType Origin; + double o2Dx, o2Dy; + + o2Dx = ((double)this->m_Size[0] - 1.) / 2.; + o2Dy = ((double)this->m_Size[1] - 1.) / 2.; + // Compute the origin (in mm) of the 2D Image + Origin[0] = -(this->m_Spacing[0]) * o2Dx; + Origin[1] = -(this->m_Spacing[1]) * o2Dy; + Origin[2] = -(this->m_SCD) ; + + return + Origin; +} + + +void DRTImageMetaInformation::SetSizeWithBounds( + double* dBounds, + SizeType MaxSize, + SpacingType Spacing) +{ + + double dXmin = dBounds[0]; + double dYmin = dBounds[1]; + double dZmin = dBounds[2]; + double dXmax = dBounds[3]; + double dYmax = dBounds[4]; + double dZmax = dBounds[5]; + + double + dNPixMin1 = 0., + dNPixMin2 = 0., + dNPixMax1 = 0., + dNPixMax2 = 0.; + + dNPixMin1 = abs(dXmin / Spacing[0]); + dNPixMin2 = abs(dYmin / Spacing[0]); + dNPixMax1 = abs(dXmax / Spacing[0]); + dNPixMax2 = abs(dYmax / Spacing[0]); + + double dPixMaxOverlap =dNPixMin1; + + if(dPixMaxOverlap > dNPixMin2){ + dPixMaxOverlap = dNPixMin2; + } + + if(dPixMaxOverlap > dNPixMax1){ + dPixMaxOverlap = dNPixMax1; + } + + if(dPixMaxOverlap > dNPixMax2){ + dPixMaxOverlap = dNPixMax2; + } + + if(floor(dPixMaxOverlap*2) > MaxSize[0]) { + m_Size[0] = MaxSize[0]; + } else { + m_Size[0] = floor(dPixMaxOverlap*2); + } + + + + dNPixMin1 = abs(dZmin / Spacing[1]); + dNPixMax1 = abs(dZmax / Spacing[1]); + + dPixMaxOverlap =dNPixMin1; + if(dPixMaxOverlap > dNPixMax1){ + dPixMaxOverlap = dNPixMax1; + } + if(floor(dPixMaxOverlap*2) > MaxSize[1]) { + m_Size[1] = MaxSize[1]; + } else { + m_Size[1] = floor(dPixMaxOverlap*2); + } + + m_Size[2] = 1; + + m_Spacing = Spacing; + + + + return; + + + +} + + DRTImageMetaInformation::PointType DRTImageMetaInformation::GetLPStoProjectionGeoLPSOffset() { @@ -137,19 +229,117 @@ DRTImageMetaInformation ::GetCOV() PointType LastVoxelPosition = Dest; - LastVoxelPosition [0] += m_Origin[0]; - LastVoxelPosition [1] += m_Origin[1]; - LastVoxelPosition [2] += m_Origin[2]; + PointType Origin = this->GetOrigin(); + LastVoxelPosition [0] += Origin[0]; + LastVoxelPosition [1] += Origin[1]; + LastVoxelPosition [2] += 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.; + image3DCOV[0] = (Origin [0] + LastVoxelPosition[0])/2.; + image3DCOV[1] = (Origin [1] + LastVoxelPosition[1])/2.; + image3DCOV[2] = (Origin [2] + LastVoxelPosition[2])/2.; return image3DCOV; } +CTVolumeImageMetaInformation::PointType +CTVolumeImageMetaInformation::ConvertIECPointToLPSPoint( + PointType IECPoint){ + + DirectionType IECtoLPS_Directions; + IECtoLPS_Directions = this->m_LPS2IECDirections.GetTranspose(); + + PointType m_LPSPoint; + m_LPSPoint = + IECtoLPS_Directions * IECPoint; + + return + m_LPSPoint; +} + +CTVolumeImageMetaInformation::PointType +CTVolumeImageMetaInformation::GetProjectionOriginLPS( + PointType ProjectionCenter) +{ + +// PointType Origin; +// Origin = this->m_OriginLPS; +// Origin = Origin - this->m_ImportOffset; + + DirectionType IECtoLPS_Directions; + IECtoLPS_Directions = this->m_LPS2IECDirections.GetTranspose(); + + + PointType m_ProjectionOriginLPS; + m_ProjectionOriginLPS = + IECtoLPS_Directions * ProjectionCenter; + + /* in longitudinal direction (Sup-Inf), we put it in the + * middle of the volume */ + m_ProjectionOriginLPS[2] = this->GetCOV()[2] - this->m_ImportOffset[2]; + + return + m_ProjectionOriginLPS; + +} + +CTVolumeImageMetaInformation::PointType +CTVolumeImageMetaInformation::GetProjectionOriginLPSZero( + PointType ProjectionCenter) +{ + + PointType Origin; + Origin = this->m_OriginLPS; + Origin = Origin - this->m_ImportOffset; + + + DirectionType IECtoLPS_Directions; + IECtoLPS_Directions = this->m_LPS2IECDirections.GetTranspose(); + + PointType m_ProjectionOriginLPS; + m_ProjectionOriginLPS = + IECtoLPS_Directions * ProjectionCenter; + /* in longitudinal direction (Sup-Inf), we put it in the + * middle of the volume */ + m_ProjectionOriginLPS[2] = this->GetCOV()[2] - this->m_ImportOffset[2]; + + + return + m_ProjectionOriginLPS - Origin; + +} + +const std::vector +CTVolumeImageMetaInformation::CalculateRegions(PointType pProjectionOriginLPS) +{ + + /* 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 = + (m_ImageDirections * Dest); + + LastVoxelPosition [0] += m_OriginLPS[0]; + LastVoxelPosition [1] += m_OriginLPS[1]; + LastVoxelPosition [2] += m_OriginLPS[2]; + + std::vector vBounds; + vBounds.clear(); + vBounds.push_back( pProjectionOriginLPS[0] - LastVoxelPosition [0] ); + vBounds.push_back( pProjectionOriginLPS[1] - LastVoxelPosition [1] ); + vBounds.push_back( pProjectionOriginLPS[2] - LastVoxelPosition [2] ); + vBounds.push_back( pProjectionOriginLPS[0] - m_OriginLPS [0] ); + vBounds.push_back( pProjectionOriginLPS[1] - m_OriginLPS [1] ); + vBounds.push_back( pProjectionOriginLPS[2] - m_OriginLPS [2] ); + + return + vBounds; + +} CTVolumeImageMetaInformation:: CTVolumeImageMetaInformation(){ @@ -250,7 +440,17 @@ DRTProjectionGeometryImageMetaInformation(){ this->m_ProjectionAngle2IEC = 0.; - this->m_SCD = 0.; + this->m_ProjectionAngle1OffsetIEC = 0.; + + this->m_ProjectionAngle2OffsetIEC = 0.; + + this->m_SCD1 = 0.; + + this->m_SCD2 = 0.; + + this->m_SCD1Offset = 0.; + + this->m_SCD2Offset = 0.; this->m_IntensityThreshold=0.; @@ -266,6 +466,10 @@ DRTProjectionGeometryImageMetaInformation(){ this->m_IECS2IECScannerR.Fill(0.); + this->m_ProjectionCenterOffset1.Fill(0.); + + this->m_ProjectionCenterOffset2.Fill(0.); + this->m_ProjectionCenter.Fill(0.); } diff --git a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h index b270e7b..5ec3556 100644 --- a/reg23Topograms/itkDTRrecon/DRTMetaInformation.h +++ b/reg23Topograms/itkDTRrecon/DRTMetaInformation.h @@ -117,8 +117,8 @@ public: itkSetMacro(Spacing,SpacingType); itkGetMacro(Spacing,SpacingType); - itkSetMacro(Origin,PointType); - itkGetMacro(Origin,PointType); + //itkSetMacro(Origin,PointType); + //itkGetMacro(Origin,PointType); itkSetMacro(OriginLPS,PointType); itkGetMacro(OriginLPS,PointType); @@ -141,6 +141,10 @@ public: PointType GetCOV(); PointType GetLPStoProjectionGeoLPSOffset(); + void SetSizeWithBounds(double*, SizeType, SpacingType); + + PointType GetOrigin(); + protected: SizeType @@ -150,7 +154,7 @@ protected: m_Spacing; PointType - m_Origin, + //m_Origin, m_OriginLPS, m_ProjectionOriginLPS, m_ProjectionOriginLPSZero; @@ -222,6 +226,18 @@ public: PointType GetCOV(); PointType GetOriginWOffset(); + const std::vector CalculateRegions(PointType pProjectionOriginLPS); + + + PointType ConvertIECPointToLPSPoint( + PointType IECPoint); + + PointType GetProjectionOriginLPS( + PointType ProjectionCenter); + + PointType GetProjectionOriginLPSZero( + PointType ProjectionCenter); + protected: tPatOrientation @@ -283,8 +299,23 @@ public: itkSetMacro(ProjectionAngle2IEC, double); itkGetMacro(ProjectionAngle2IEC, double); - itkSetMacro(SCD, double); - itkGetMacro(SCD, double); + itkSetMacro(ProjectionAngle1OffsetIEC, double); + itkGetMacro(ProjectionAngle1OffsetIEC, double); + + itkSetMacro(ProjectionAngle2OffsetIEC, double); + itkGetMacro(ProjectionAngle2OffsetIEC, double); + + itkSetMacro(SCD1, double); + itkGetMacro(SCD1, double); + + itkSetMacro(SCD2, double); + itkGetMacro(SCD2, double); + + itkSetMacro(SCD1Offset, double); + itkGetMacro(SCD1Offset, double); + + itkSetMacro(SCD2Offset, double); + itkGetMacro(SCD2Offset, double); itkSetMacro(IntensityThreshold, double); itkGetMacro(IntensityThreshold, double); @@ -307,6 +338,12 @@ public: itkSetMacro(IECS2IECScannerR, PointType); itkGetMacro(IECS2IECScannerR, PointType); + itkSetMacro(ProjectionCenterOffset1, PointType); + itkGetMacro(ProjectionCenterOffset1, PointType); + + itkSetMacro(ProjectionCenterOffset2, PointType); + itkGetMacro(ProjectionCenterOffset2, PointType); + itkSetMacro(ProjectionCenter, PointType); itkGetMacro(ProjectionCenter, PointType); @@ -315,7 +352,12 @@ protected: double m_ProjectionAngle1IEC, m_ProjectionAngle2IEC, - m_SCD, + m_ProjectionAngle1OffsetIEC, + m_ProjectionAngle2OffsetIEC, + m_SCD1, + m_SCD2, + m_SCD1Offset, + m_SCD2Offset, m_IntensityThreshold; @@ -337,7 +379,9 @@ protected: PointType /* center of projection in an IEC reference at * Patient Origin of fixed images. Positioning scanner */ - m_ProjectionCenter; + m_ProjectionCenter, + m_ProjectionCenterOffset1, + m_ProjectionCenterOffset2; /** Default Constructor **/ diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index d8b18c2..59ac4d9 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -85,8 +85,10 @@ itkImageProcessor::itkImageProcessor() { // std::cout<<"itkImageProcessor::itkImageProcessor contructor"<SetIdentity(); + transform1 = TransformType::New(); + transform1->SetIdentity(); + transform2 = TransformType::New(); + transform2->SetIdentity(); interpolator1 = InterpolatorType::New(); interpolator2 = InterpolatorType::New(); @@ -143,7 +145,8 @@ itkImageProcessor::itkImageProcessor() /* Initialise the projection geoemtry with defaults */ m_DRTGeometryMetaInfo = DRTProjectionGeometryImageMetaInformation::New(); - m_DRTGeometryMetaInfo->SetSCD(570.); + m_DRTGeometryMetaInfo->SetSCD1(570.); + m_DRTGeometryMetaInfo->SetSCD2(570.); m_DRTGeometryMetaInfo->SetProjectionAngle1IEC(180.); m_DRTGeometryMetaInfo->SetProjectionAngle2IEC(90.); m_DRTGeometryMetaInfo->SetIntensityThreshold(-300.); @@ -153,6 +156,7 @@ itkImageProcessor::itkImageProcessor() Point3D[1]=0.; Point3D[2]=175.; m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D); + m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D); ImageType3D::SizeType ImageSize; ImageSize[0]=512; @@ -172,7 +176,8 @@ itkImageProcessor::itkImageProcessor() Point3D[2]=0.; m_DRTGeometryMetaInfo->SetIECS2IECScannerT(Point3D); m_DRTGeometryMetaInfo->SetIECS2IECScannerR(Point3D); - + m_DRTGeometryMetaInfo->SetProjectionCenterOffset1(Point3D); + m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(Point3D); // // TEST MAPPING @@ -229,14 +234,25 @@ void itkImageProcessor::SetIntensityThreshold(double dT) m_DRTGeometryMetaInfo->SetIntensityThreshold(dT); } -void itkImageProcessor::SetSCD(double dDist) +void itkImageProcessor::SetSCDOffset(double dOff1, double dOff2) { - m_DRTGeometryMetaInfo->SetSCD(dDist); + m_DRTGeometryMetaInfo->SetSCD1Offset(dOff1); + m_DRTGeometryMetaInfo->SetSCD2Offset(dOff2); } -double itkImageProcessor::GetSCD(){ +void itkImageProcessor::SetSCD(double dDist1, double dDist2) +{ + m_DRTGeometryMetaInfo->SetSCD1(dDist1); + m_DRTGeometryMetaInfo->SetSCD2(dDist2); +} + +double itkImageProcessor::GetSCD1(){ return - m_DRTGeometryMetaInfo->GetSCD(); + m_DRTImage1MetaInfo->GetSCD(); +} +double itkImageProcessor::GetSCD2(){ + return + m_DRTImage1MetaInfo ->GetSCD(); } void itkImageProcessor::SetCustom_ImportTransform(double dTx, @@ -247,6 +263,10 @@ void itkImageProcessor::SetCustom_ImportTransform(double dTx, double dRz) { + if(m_DRTGeometryMetaInfo == NULL){ + //todo + } + ImageType3D::PointType Punto; Punto[0] = dTx; @@ -263,17 +283,45 @@ void itkImageProcessor::SetCustom_ImportTransform(double dTx, } -void itkImageProcessor::SetCustom_ProjectionCenterFixedAxes_IEC(double dX, - double dY, - double dZ) +void itkImageProcessor::SetCustom_ProjectionCenterOffsetFixedAxes_IEC( + double dX1, + double dY1, + double dZ1, + double dX2, + double dY2, + double dZ2) { typedef itk::Point PointType; PointType m_point; - m_point [0] = dX; - m_point [1] = dY; - m_point [2] = dZ; + m_point [0] = dX1; + m_point [1] = dY1; + m_point [2] = dZ1; + + m_DRTGeometryMetaInfo->SetProjectionCenterOffset1(m_point); + + m_point [0] = dX2; + m_point [1] = dY2; + m_point [2] = dZ2; + + m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(m_point); + +} + + +void itkImageProcessor::SetCustom_ProjectionCenterFixedAxes_IEC( + double dX1, + double dY1, + double dZ1) +{ + + typedef itk::Point PointType; + PointType m_point; + + m_point [0] = dX1; + m_point [1] = dY1; + m_point [2] = dZ1; m_DRTGeometryMetaInfo->SetProjectionCenter(m_point); @@ -480,36 +528,56 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) PointOffset.Fill(0.); m_CTMetaInfo->SetImportOffset(PointOffset); + + /* initialise DRT meta */ m_DRTImage1MetaInfo = DRTImageMetaInformation::New(); m_DRTImage1MetaInfo->SetProjectionAngleLPS( this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), - m_DRTGeometryMetaInfo->GetProjectionAngle1IEC()) + m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() + + m_DRTGeometryMetaInfo->GetProjectionAngle1OffsetIEC()) ); + m_DRTImage1MetaInfo->SetSCD( - m_DRTGeometryMetaInfo->GetSCD()); + m_DRTGeometryMetaInfo->GetSCD1() + + m_DRTGeometryMetaInfo->GetSCD1Offset() + ); /* Calculate projection angle IEC to LPS */ - m_DRTImage2MetaInfo = DRTImageMetaInformation::New(); m_DRTImage2MetaInfo->SetProjectionAngleLPS( this->CalcProjectionAngleLPS( m_CTMetaInfo->GetPatientOrientation(), - m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) + m_DRTGeometryMetaInfo->GetProjectionAngle2IEC() + + m_DRTGeometryMetaInfo->GetProjectionAngle2OffsetIEC()) ); + m_DRTImage2MetaInfo->SetSCD( - m_DRTGeometryMetaInfo->GetSCD()); + m_DRTGeometryMetaInfo->GetSCD2()+ + m_DRTGeometryMetaInfo->GetSCD2Offset() + ); + std::cout<< "///////////////// 3D VOLUME BEG ///////////////" <GetOriginLPS() <GetCOV() <GetImportOffset() <GetLPS2IECDirections() <GetSize() <GetSpacing() <GetProjectionAngleLPS() <GetSCD() <GetProjectionAngleLPS() <GetSCD() <UpdateProjectionGeometryMeta(); - // To simply Siddon-Jacob's fast ray-tracing algorithm, we force the origin of the CT image // 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. @@ -534,64 +602,6 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) } -itkImageProcessor::ImageType3D::PointType -itkImageProcessor::CalculateDRTOrigin( - ImageType3D::SizeType m_ImageSize, - ImageType3D::SpacingType m_ImageResolution, - double dSCD - ){ - - 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 Origin; - Origin = m_VolumeOriginLPS; - Origin = Origin - m_ImportOffset; - - - InternalImageType::DirectionType IECtoLPS_Directions; - IECtoLPS_Directions = m_VolumeLPS2IEC.GetTranspose(); - - 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) - { - return - m_ProjectionOriginLPS - Origin; - } else { - return - m_ProjectionOriginLPS; - } -} - - const std::vector itkImageProcessor::GetRTImportOffset() { std::vector vOffset; @@ -695,6 +705,12 @@ itkImageProcessor::CalcProjectionAngleLPS( } +void itkImageProcessor::SetProjectionAngleOffsetIEC(double dOff1, double dOff2) +{ + m_DRTGeometryMetaInfo->SetProjectionAngle1OffsetIEC(dOff1); + m_DRTGeometryMetaInfo->SetProjectionAngle2OffsetIEC(dOff2); +} + void itkImageProcessor::SetProjectionAngle1IEC(double dGantryAngle) { @@ -922,6 +938,15 @@ int itkImageProcessor::load2D(const char * pcFName){ m_Duplicator->Update(); + std::cout<< " ////////////// 2D Topo BEG " <GetOutput()->GetDirection()<GetOutput()->GetOrigin()<GetOutput()->GetSpacing()<GetOutput()->GetLargestPossibleRegion().GetSize()<GetCenter() "<< m_UserTransform->GetCenter() <GetTranslation() "<< m_UserTransform->GetTranslation() <GetAngleX() "<< m_UserTransform->GetAngleX() <SetIdentity(); m_OutputTransform->SetTranslation( - //m_OffsetTransform->GetTranslation() + + m_OffsetTransform->GetTranslation() + m_UserTransform->GetTranslation() ); m_OutputTransform->SetRotation( - //m_OffsetTransform->GetAngleX() + + m_OffsetTransform->GetAngleX() + m_UserTransform->GetAngleX(), - //m_OffsetTransform->GetAngleY() + + m_OffsetTransform->GetAngleY() + m_UserTransform->GetAngleY(), - //m_OffsetTransform->GetAngleZ() + + m_OffsetTransform->GetAngleZ() + m_UserTransform->GetAngleZ() ); @@ -1138,93 +1155,184 @@ void itkImageProcessor::InitializeProjector() //TODO } + const double dtr = (atan(1.0) * 4.0) / 180.0; + + ImageType3D::PointType ZeroPoint; ZeroPoint.Fill(0.); - std::cout<< "itkImageProcessor::InitializeProjector()" <GetLPS2IECDirections().GetTranspose(); - TransformType::Pointer CurrTransform = - CalculateInternalTransform( - ZeroPoint , - ZeroPoint , - m_TransformMetaInfo->GetTranslations(), - m_TransformMetaInfo->GetRotations(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), - ZeroPoint, - IECtoLPS_Directions - ); + TransformType::Pointer CurrTransform; + //CurrTransform->SetComputeZYX(true); + //CurrTransform->SetIdentity(); - transform->SetComputeZYX(true); - transform->SetIdentity(); + if(m_RTMetaInfo == NULL) + { - transform->SetTranslation( + //AAAAAA TODOOOOO CERCA + ImageType3D::PointType pFakeIsoLPS = + m_DRTImage1MetaInfo->GetProjectionOriginLPS(); + + //ImageType3D::PointType pFakeIsoLPS = + // m_CTMetaInfo->GetProjectionOriginLPS( + // m_DRTGeometryMetaInfo->GetProjectionCenter1()); + + + CurrTransform = + CalculateInternalTransform( + ZeroPoint, + ZeroPoint, + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage1MetaInfo->GetProjectionOriginLPS(), + pFakeIsoLPS, + ZeroPoint, + IECtoLPS_Directions + ); + } else { + + CurrTransform = + CalculateInternalTransform( + ZeroPoint , + m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage1MetaInfo->GetProjectionOriginLPS(), + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), + m_CTMetaInfo->GetImportOffset(), + IECtoLPS_Directions + ); + + } + + + transform1->SetComputeZYX(true); + transform1->SetIdentity(); + + transform1->SetTranslation( CurrTransform->GetTranslation()); - transform->SetRotation( + transform1->SetRotation( CurrTransform->GetAngleX(), CurrTransform->GetAngleY(), CurrTransform->GetAngleZ() ); -// transform->SetCenter( -// CurrTransform->GetCenter() -// ); - std::cout<< "itkImageProcessor::InitializeProjector()" <GetTranslations()[0]; //TZero[0]; -// translation[1] = m_TransformMetaInfo->GetTranslations()[1]; //TZero[1]; -// translation[2] = m_TransformMetaInfo->GetTranslations()[2]; //TZero[2]; -// transform->SetTranslation(translation); - const double dtr = (atan(1.0) * 4.0) / 180.0; -// transform->SetRotation( -// dtr * m_TransformMetaInfo->GetRotations()[0],// RZero[0], -// dtr * m_TransformMetaInfo->GetRotations()[1], // RZero[1], -// dtr * m_TransformMetaInfo->GetRotations()[2]); //RZero[2]); - -// if(m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() != -// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ) { -// // TODO -// } - - transform->SetCenter( + transform1->SetCenter( m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() ); - transform->Print(std::cout); + transform1->Print(std::cout); // 2D Image 1 interpolator1->SetProjectionAngle( - dtr * - m_DRTImage1MetaInfo->GetProjectionAngleLPS() ); + dtr * m_DRTImage1MetaInfo->GetProjectionAngleLPS() ); interpolator1->SetFocalPointToIsocenterDistance( m_DRTImage1MetaInfo->GetSCD()); interpolator1->SetThreshold( - m_DRTGeometryMetaInfo->GetIntensityThreshold() - ); - interpolator1->SetTransform(transform); + m_DRTGeometryMetaInfo->GetIntensityThreshold() ); + interpolator1->SetTransform(transform1); interpolator1->Initialize(); + + + + + + + //CurrTransform->SetComputeZYX(true); + //CurrTransform->SetIdentity(); + + if(m_RTMetaInfo == NULL) + { + ImageType3D::PointType pFakeIsoLPS = + m_DRTImage2MetaInfo->GetProjectionOriginLPS(); + + //ImageType3D::PointType pFakeIsoLPS = + // m_CTMetaInfo->GetProjectionOriginLPS( + // m_DRTGeometryMetaInfo->GetProjectionCenter2()); + + + + CurrTransform = + CalculateInternalTransform( + ZeroPoint, + ZeroPoint, + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + pFakeIsoLPS, + ZeroPoint, + IECtoLPS_Directions + ); + } else { + + CurrTransform = + CalculateInternalTransform( + ZeroPoint , + m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), + m_CTMetaInfo->GetImportOffset(), + IECtoLPS_Directions + ); + + } + + + transform2->SetComputeZYX(true); + transform2->SetIdentity(); + + transform2->SetTranslation( + CurrTransform->GetTranslation()); + transform2->SetRotation( + CurrTransform->GetAngleX(), + CurrTransform->GetAngleY(), + CurrTransform->GetAngleZ() + ); + + + transform2->SetCenter( + m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ); + //transform2->Print(std::cout); + // 2D Image 2 interpolator2->SetProjectionAngle( - dtr * - m_DRTImage2MetaInfo->GetProjectionAngleLPS() ); + dtr * m_DRTImage2MetaInfo->GetProjectionAngleLPS() ); interpolator2->SetFocalPointToIsocenterDistance( - m_DRTImage2MetaInfo->GetSCD()); + m_DRTImage2MetaInfo->GetSCD() ); interpolator2->SetThreshold( - m_DRTGeometryMetaInfo->GetIntensityThreshold() - ); - interpolator2->SetTransform(transform); + m_DRTGeometryMetaInfo->GetIntensityThreshold() ); + interpolator2->SetTransform(transform2); interpolator2->Initialize(); + + +// // 2D Image 2 +// interpolator2->SetProjectionAngle( +// dtr * +// m_DRTImage2MetaInfo->GetProjectionAngleLPS() ); +// interpolator2->SetFocalPointToIsocenterDistance( +// m_DRTImage2MetaInfo->GetSCD()); +// interpolator2->SetThreshold( +// m_DRTGeometryMetaInfo->GetIntensityThreshold() +// ); +// interpolator2->SetTransform(transform); +// interpolator2->Initialize(); + + + + + resampleFilter1 = ResampleFilterType::New(); resampleFilter1->SetInput( image3DIn); @@ -1235,7 +1343,7 @@ void itkImageProcessor::InitializeProjector() // The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance // have been set before registration. Here we only need to replace the initial // transform with the final transform. - interpolator1->SetTransform(transform); + interpolator1->SetTransform(transform1); interpolator1->Initialize(); resampleFilter1->SetInterpolator(interpolator1); @@ -1255,7 +1363,7 @@ void itkImageProcessor::InitializeProjector() // The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance // have been set before registration. Here we only need to replace the initial // transform with the final transform. - interpolator2->SetTransform(transform); + interpolator2->SetTransform(transform2); interpolator2->Initialize(); resampleFilter2->SetInterpolator(interpolator2); @@ -1296,7 +1404,7 @@ void itkImageProcessor::InitializeProjector() filter1->Update(); filter2->Update(); - std::cout<< "itkImageProcessor::InitializeProjector() AFTER UPDATE" <GetOutput(); imageDRT2In = filter2->GetOutput(); @@ -1346,7 +1454,7 @@ void itkImageProcessor::loadRTPlanInfo( //TODO } - //m_RTMetaInfo = RTGeometryMetaInformation::New(); + m_RTMetaInfo = RTGeometryMetaInformation::New(); ImageType3D::PointType Point; @@ -1367,7 +1475,8 @@ void itkImageProcessor::loadRTPlanInfo( m_RTMetaInfo->GetIsocenterIECS(), m_CTMetaInfo->GetLPS2IECDirections(), m_RTMetaInfo->GetIsocenterLPS(), - m_DRTGeometryMetaInfo->GetIECS2IECScannerT() + m_DRTGeometryMetaInfo->GetIECS2IECScannerT(), + m_DRTGeometryMetaInfo->GetIECS2IECScannerR() ) ); @@ -1375,7 +1484,7 @@ void itkImageProcessor::loadRTPlanInfo( <GetImportOffset() <UpdateProjectionGeometryMeta(); - + this->InitializeProjector(); } @@ -1397,104 +1506,209 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ //TODO. } + + // FIRST + + ImageType3D::PointType NominalIsocenter_wZcorrectionLPS; + + NominalIsocenter_wZcorrectionLPS = + m_CTMetaInfo->GetProjectionOriginLPS( + m_DRTGeometryMetaInfo->GetProjectionCenter() + ); + + ImageType3D::PointType NominalIsocenterZero_wZcorrectionLPS; + + NominalIsocenterZero_wZcorrectionLPS = + m_CTMetaInfo->GetProjectionOriginLPSZero( + m_DRTGeometryMetaInfo->GetProjectionCenter() + ); + + std::cout<< "///////////////// PGEOM META BEG ///////////////" <GetProjectionCenter() <ConvertIECPointToLPSPoint( + m_DRTGeometryMetaInfo->GetProjectionCenterOffset1()); + + + + ImageType3D::PointType CalibratedIsocenterLPS; + + CalibratedIsocenterLPS[0] = NominalIsocenter_wZcorrectionLPS[0] + + IsocenterOffsetLPS[0]; + CalibratedIsocenterLPS[1] = NominalIsocenter_wZcorrectionLPS[1] + + IsocenterOffsetLPS[1]; + CalibratedIsocenterLPS[2] = NominalIsocenter_wZcorrectionLPS[2] + + IsocenterOffsetLPS[2]; + + + ImageType3D::PointType CalibratedIsocenterZeroLPS; + + CalibratedIsocenterZeroLPS[0] = NominalIsocenterZero_wZcorrectionLPS[0] + + IsocenterOffsetLPS[0]; + CalibratedIsocenterZeroLPS[1] = NominalIsocenterZero_wZcorrectionLPS[1] + + IsocenterOffsetLPS[1]; + CalibratedIsocenterZeroLPS[2] = NominalIsocenterZero_wZcorrectionLPS[2] + + IsocenterOffsetLPS[2]; + + + + m_DRTImage1MetaInfo->SetProjectionOriginLPS( - CalculateProjectionCenterLPS( - m_CTMetaInfo->GetImportOffset(), - m_CTMetaInfo->GetLPS2IECDirections(), - m_CTMetaInfo->GetCOV(), - m_CTMetaInfo->GetOriginLPS(), - m_DRTGeometryMetaInfo->GetProjectionCenter(), - false - ) - ); + CalibratedIsocenterLPS); + m_DRTImage1MetaInfo->SetProjectionOriginLPSZero( - CalculateProjectionCenterLPS( - m_CTMetaInfo->GetImportOffset(), - m_CTMetaInfo->GetLPS2IECDirections(), - m_CTMetaInfo->GetCOV(), - m_CTMetaInfo->GetOriginLPS(), - m_DRTGeometryMetaInfo->GetProjectionCenter(), - true - ) + CalibratedIsocenterZeroLPS); + + // This is based on the calibrated iso to be sure... + std::vector vBounds = m_CTMetaInfo->CalculateRegions( + NominalIsocenter_wZcorrectionLPS//m_DRTImage1MetaInfo->GetProjectionOriginLPS() ); - //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() - ) + + std::cout<<"Bounds1: " + <SetSizeWithBounds( + vBounds.data(), + m_DRTGeometryMetaInfo->GetDRT1Size(), + m_DRTGeometryMetaInfo->GetDRT1Spacing() ); + + + // This HAS TO be calculated from the nominal isocenter. m_DRTImage1MetaInfo->SetOriginLPS( CalcDRTImageOrigin( m_DRTImage1MetaInfo->GetOrigin(), m_DRTImage1MetaInfo->GetCOV(), - m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - m_DRTImage1MetaInfo->GetProjectionAngleLPS(), + NominalIsocenter_wZcorrectionLPS,//m_DRTImage1MetaInfo->GetProjectionOriginLPS(), + this->CalcProjectionAngleLPS( + m_CTMetaInfo->GetPatientOrientation(), + m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() ), + //m_DRTGeometryMetaInfo->GetProjectionAngle1IEC(),// m_DRTImage1MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS ) ); + + + m_DRTImage1MetaInfo->SetImageDirectionsLPS( CalcDRTImageDirections( - m_DRTImage1MetaInfo->GetProjectionAngleLPS(), + this->CalcProjectionAngleLPS( + m_CTMetaInfo->GetPatientOrientation(), + m_DRTGeometryMetaInfo->GetProjectionAngle1IEC()), + //m_DRTGeometryMetaInfo->GetProjectionAngle1IEC(),//m_DRTImage1MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS) ); + std::cout<<"IsocenterOffsetLPS"<< IsocenterOffsetLPS <GetSize()"<< m_DRTImage1MetaInfo->GetSize() <GetSpacing()"<< m_DRTImage1MetaInfo->GetSpacing() <GetOrigin()"<< m_DRTImage1MetaInfo->GetOrigin() <GetCOV()"<< m_DRTImage1MetaInfo->GetCOV() <GetOriginLPS()"<< m_DRTImage1MetaInfo->GetOriginLPS() <GetProjectionAngleLPS()"<GetProjectionAngleLPS()<GetImageDirectionsLPS()"<< m_DRTImage1MetaInfo->GetImageDirectionsLPS() <ConvertIECPointToLPSPoint( + m_DRTGeometryMetaInfo->GetProjectionCenterOffset2()); + + CalibratedIsocenterLPS[0] = NominalIsocenter_wZcorrectionLPS[0] + + IsocenterOffsetLPS[0]; + CalibratedIsocenterLPS[1] = NominalIsocenter_wZcorrectionLPS[1] + + IsocenterOffsetLPS[1]; + CalibratedIsocenterLPS[2] = NominalIsocenter_wZcorrectionLPS[2] + + IsocenterOffsetLPS[2]; + + + CalibratedIsocenterZeroLPS[0] = NominalIsocenterZero_wZcorrectionLPS[0] + + IsocenterOffsetLPS[0]; + CalibratedIsocenterZeroLPS[1] = NominalIsocenterZero_wZcorrectionLPS[1] + + IsocenterOffsetLPS[1]; + CalibratedIsocenterZeroLPS[2] = NominalIsocenterZero_wZcorrectionLPS[2] + + IsocenterOffsetLPS[2]; + + + m_DRTImage2MetaInfo->SetProjectionOriginLPS( - CalculateProjectionCenterLPS( - m_CTMetaInfo->GetImportOffset(), - m_CTMetaInfo->GetLPS2IECDirections(), - m_CTMetaInfo->GetCOV(), - m_CTMetaInfo->GetOriginLPS(), - m_DRTGeometryMetaInfo->GetProjectionCenter(), - false - ) - ); + CalibratedIsocenterLPS); +// m_CTMetaInfo->GetProjectionOriginLPS( +// m_DRTGeometryMetaInfo->GetProjectionCenter2()) +// ); + m_DRTImage2MetaInfo->SetProjectionOriginLPSZero( - CalculateProjectionCenterLPS( - m_CTMetaInfo->GetImportOffset(), - m_CTMetaInfo->GetLPS2IECDirections(), - m_CTMetaInfo->GetCOV(), - m_CTMetaInfo->GetOriginLPS(), - m_DRTGeometryMetaInfo->GetProjectionCenter(), - true - ) + CalibratedIsocenterZeroLPS); +// m_CTMetaInfo->GetProjectionOriginLPSZero( +// m_DRTGeometryMetaInfo->GetProjectionCenter2()) +// ); + + vBounds.clear(); + vBounds = m_CTMetaInfo->CalculateRegions( + NominalIsocenter_wZcorrectionLPS//m_DRTImage2MetaInfo->GetProjectionOriginLPS() ); - //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() - ) + + std::cout<<"Bounds2: " + <SetSizeWithBounds( + vBounds.data(), + m_DRTGeometryMetaInfo->GetDRT2Size(), + m_DRTGeometryMetaInfo->GetDRT2Spacing() ); + m_DRTImage2MetaInfo->SetOriginLPS( CalcDRTImageOrigin( m_DRTImage2MetaInfo->GetOrigin(), m_DRTImage2MetaInfo->GetCOV(), - m_DRTImage2MetaInfo->GetProjectionOriginLPS(), - m_DRTImage2MetaInfo->GetProjectionAngleLPS(), + NominalIsocenter_wZcorrectionLPS,//m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + this->CalcProjectionAngleLPS( + m_CTMetaInfo->GetPatientOrientation(), + m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , + //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS ) ); + m_DRTImage2MetaInfo->SetImageDirectionsLPS( CalcDRTImageDirections( - m_DRTImage2MetaInfo->GetProjectionAngleLPS(), + this->CalcProjectionAngleLPS( + m_CTMetaInfo->GetPatientOrientation(), + m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , + //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS) ); - std::cout<< "Projection Origin LPS" <GetProjectionOriginLPS() <GetOriginLPS() <GetImageDirectionsLPS() <GetProjectionOriginLPS() <GetOriginLPS() <GetImageDirectionsLPS() <GetSize()"<< m_DRTImage2MetaInfo->GetSize() <GetSpacing()"<< m_DRTImage2MetaInfo->GetSpacing() <GetOrigin()"<< m_DRTImage2MetaInfo->GetOrigin() <GetCOV()"<< m_DRTImage2MetaInfo->GetCOV() <GetOriginLPS()"<< m_DRTImage2MetaInfo->GetOriginLPS() <GetProjectionAngleLPS()"<GetProjectionAngleLPS()<GetImageDirectionsLPS()"<< m_DRTImage2MetaInfo->GetImageDirectionsLPS() <SetComputeZYX(true); - transform->SetIdentity(); + transform1->SetComputeZYX(true); + transform1->SetIdentity(); - InternalImageType::DirectionType IECtoLPS_Directions; + InternalImageType::DirectionType IECtoLPS_Directions; - IECtoLPS_Directions = - m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); + IECtoLPS_Directions = + m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); + TransformType::Pointer CurrTransform; + //CurrTransform->SetComputeZYX(true); + //CurrTransform->SetIdentity(); - TransformType::Pointer CurrTransform = + if(m_RTMetaInfo == NULL) + { + + ImageType3D::PointType pFakeIsoLPS = + m_DRTImage1MetaInfo->GetProjectionOriginLPS(); + //m_CTMetaInfo->GetProjectionOriginLPS( + // m_DRTGeometryMetaInfo->GetProjectionCenter()); + + std::cout<<"pFakeIsoLPS: "<GetTranslations(), m_TransformMetaInfo->GetRotations(), m_DRTImage1MetaInfo->GetProjectionOriginLPS(), - m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), + pFakeIsoLPS, ZeroPoint, IECtoLPS_Directions ); + } else { - transform->SetComputeZYX(true); - transform->SetIdentity(); + CurrTransform = + CalculateInternalTransform( + ZeroPoint , + m_DRTGeometryMetaInfo->GetIECS2IECScannerR(), + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage1MetaInfo->GetProjectionOriginLPS(), + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), + m_CTMetaInfo->GetImportOffset(), + IECtoLPS_Directions + ); - transform->SetTranslation( + } + + + + transform1->SetComputeZYX(true); + transform1->SetIdentity(); + + transform1->SetTranslation( + CurrTransform->GetTranslation()); + transform1->SetRotation( + CurrTransform->GetAngleX(), + CurrTransform->GetAngleY(), + CurrTransform->GetAngleZ() + ); + + std::cout<< "itkImageProcessor::GetProjectionImages" <SetCenter( + m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()); + + std::cout<<"Angle1 "<GetProjectionAngleLPS() <GetProjectionOriginLPSZero() <GetProjectionOriginLPS()<GetProjectionAngleLPS() <GetProjectionOriginLPSZero() <GetProjectionOriginLPS()<Print(std::cout); + // The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance + // have been set before registration. Here we only need to replace the initial + // transform with the final transform. + interpolator1->SetTransform(transform1); + interpolator1->Initialize(); + + + + + + transform2->SetComputeZYX(true); + transform2->SetIdentity(); + + //CurrTransform->SetComputeZYX(true); + //CurrTransform->SetIdentity(); + + if(m_RTMetaInfo == NULL) + { + + ImageType3D::PointType pFakeIsoLPS = + m_DRTImage2MetaInfo->GetProjectionOriginLPS(); + //m_CTMetaInfo->GetProjectionOriginLPS( + // m_DRTGeometryMetaInfo->GetProjectionCenter()); + + std::cout<<"pFakeIsoLPS: "<GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + pFakeIsoLPS, + ZeroPoint, + IECtoLPS_Directions + ); + + std::cout<<"pFakeIsoLPS: "<GetIECS2IECScannerR(), + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetImportOffset(), + m_CTMetaInfo->GetImportOffset(), + IECtoLPS_Directions + ); + + } + + + + transform2->SetComputeZYX(true); + transform2->SetIdentity(); + + transform2->SetTranslation( CurrTransform->GetTranslation()); - transform->SetRotation( + transform2->SetRotation( CurrTransform->GetAngleX(), CurrTransform->GetAngleY(), CurrTransform->GetAngleZ() ); -// transform->SetCenter( -// CurrTransform->GetCenter() -// ); - std::cout<< "itkImageProcessor::GetProjectionImages" <SetComputeZYX(true); -// transform->SetIdentity(); -// TransformType::OutputVectorType translation; -// translation[0] = m_TransformMetaInfo->GetTranslations()[0]; //TZero[0]; -// translation[1] = m_TransformMetaInfo->GetTranslations()[1]; //TZero[1]; -// translation[2] = m_TransformMetaInfo->GetTranslations()[2]; //TZero[2]; - -// transform->SetTranslation(translation); -// transform->SetRotation( -// dtr * m_TransformMetaInfo->GetRotations()[0],// RZero[0], -// dtr * m_TransformMetaInfo->GetRotations()[1], // RZero[1], -// dtr * m_TransformMetaInfo->GetRotations()[2]); //RZero[2]); -// if(m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() != -// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ) { -// // TODO -// } + transform2 ->SetCenter( + m_DRTImage2MetaInfo->GetProjectionOriginLPSZero()); + + + //transform2 ->Print(std::cout); - transform ->SetCenter( - m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()); - //transform ->Print(std::cout); - // The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance - // have been set before registration. Here we only need to replace the initial - // transform with the final transform. - interpolator1->SetTransform(transform); - interpolator1->Initialize(); // // The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance // // have been set before registration. Here we only need to replace the initial // // transform with the final transform. - interpolator2->SetTransform(transform); //finalTransform); + interpolator2->SetTransform(transform2); //finalTransform); interpolator2->Initialize(); + std::cout<<"pFakeIsoLPS: "<Update(); filter2->Update(); imageDRT1In = filter1->GetOutput(); imageDRT2In = filter2->GetOutput(); - std::cout<< "itkImageProcessor::GetProjectionImages" <GetInverseTransform()->GetMatrix()<GetInverseTransform()->GetTranslation()<SetComputeZYX(true); + InputTransform->SetIdentity(); + + TransformType::OutputVectorType translation; + translation[0] = IEC2DCMMapT[0]; + translation[1] = IEC2DCMMapT[1]; + translation[2] = IEC2DCMMapT[2]; + + InputTransform->SetTranslation(translation); + + const double dtr = (atan(1.0) * 4.0) / 180.0; + InputTransform->SetRotation( + dtr * IEC2DCMMapR[0], + dtr * IEC2DCMMapR[1], + dtr * IEC2DCMMapR[2]); + + ImageType3D::PointType m_TransformOrigin; + m_TransformOrigin.Fill(0.); + InputTransform->SetCenter( + m_TransformOrigin ); + + ImageType3D::PointType IsoSupport_IECPos = + InputTransform->TransformPoint(rtCouchOffset ); + InternalImageType::DirectionType VolumeIECtoLPS; VolumeIECtoLPS = VolumeLPStoIEC.GetTranspose(); - ImageType3D::PointType IsoSupport_IECPos; - IsoSupport_IECPos[0] = rtCouchOffset[0] + IEC2DCMMapT[0]; - IsoSupport_IECPos[1] = rtCouchOffset[1] + IEC2DCMMapT[1]; - IsoSupport_IECPos[2] = rtCouchOffset[2] + IEC2DCMMapT[2]; - ImageType3D::PointType IsoSupport_LPSPos = VolumeIECtoLPS * IsoSupport_IECPos; @@ -2070,26 +2406,6 @@ void itkImageProcessor::SetInitialTranslations(double dX, double dY, double dZ) //todo } - // /* input is in IEC */ - // ImageType3D::PointType IECTranslations; - // IECTranslations[0] = dX; - // IECTranslations[1] = dY; - // IECTranslations[2] = dZ; - - // InternalImageType::DirectionType IECtoLPS_Directions; - - // IECtoLPS_Directions = - // m_CTMetaInfo->GetLPS2IECDirections().GetTranspose(); - - // ImageType3D::PointType LPSTranslations = - // IECtoLPS_Directions * IECTranslations; - - // // TZero[0]= LPSTranslations[0]; - // // TZero[1]= LPSTranslations[1]; - // // TZero[2]= LPSTranslations[2]; - - // m_TransformMetaInfo->SetTranslations(LPSTranslations); - ImageType3D::PointType Translations; Translations[0]= dX; @@ -2108,10 +2424,6 @@ void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ) m_TransformMetaInfo->SetRotations(Rotations); - // RZero[0]=dX; - // RZero[1]=dY; - // RZero[2]=dZ; - } diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.h b/reg23Topograms/itkDTRrecon/itkImageProcessor.h index ee1dfa0..30a5f68 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.h +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.h @@ -80,6 +80,7 @@ public: /** Projection angles - Gantry angle IEC */ void SetProjectionAngle1IEC(double); void SetProjectionAngle2IEC(double); + void SetProjectionAngleOffsetIEC(double,double); /** Get projection angles - Gantry angle LPS * Only meaningful after a 3D Volume is loaded */ double GetProjectionAngle1LPS(); @@ -87,8 +88,15 @@ public: /** Source to axis distance - SAD * Often called source to object (SOD), (0018,1111) */ - void SetSCD(double); - double GetSCD(); + void SetSCD(double,double); + void SetSCDOffset(double,double); + + double GetSCD1(); + double GetSCD2(); + + void SetCustom_ProjectionCenterOffsetFixedAxes_IEC(double,double,double,double,double,double); + + void SetCustom_ProjectionCenterFixedAxes_IEC(double,double,double); /** Intensity threshold used for image projection, * any voxel below threshold is ignored */ @@ -99,7 +107,7 @@ public: void SetCustom_2Dsize(int,int,int,int); void SetCustom_ImportTransform(double, double, double, double, double, double); - void SetCustom_ProjectionCenterFixedAxes_IEC(double,double,double); + /** Set transform parameters for 3D Volume */ void SetInitialTranslations(double, double, double); @@ -190,22 +198,6 @@ private: using ITKtoVTKFilterType = itk::ImageToVTKImageFilter; - 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); - TransformType::Pointer CalculateInternalTransform( ImageType3D::PointType m_TranslationOffset, @@ -220,7 +212,8 @@ private: TransformType::Pointer - transform; + transform1, + transform2; InterpolatorType::Pointer interpolator1, @@ -287,7 +280,8 @@ private: ImageType3D::PointType rtCouchOffset, InternalImageType::DirectionType VolumeLPStoIEC, ImageType3D::PointType rtIsocenterLPS, - ImageType3D::PointType IEC2DCMMapT); + ImageType3D::PointType IEC2DCMMapT, + ImageType3D::PointType IEC2DCMMapR); // The ResampleImageFilter is the driving force for the projection image generation. @@ -299,9 +293,6 @@ private: filter1, filter2; - //double TZero[3], RZero[3]; - - InternalImageType::DirectionType Std_DRT2LPS;