diff --git a/reg23Topograms/itkDTRrecon/CMakeLists.txt b/reg23Topograms/itkDTRrecon/CMakeLists.txt index 147fb1e..c93893c 100644 --- a/reg23Topograms/itkDTRrecon/CMakeLists.txt +++ b/reg23Topograms/itkDTRrecon/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 2.8.8) +cmake_minimum_required(VERSION 2.8.12) SET(LIB_NAME itkDTRrecon) @@ -10,7 +10,6 @@ find_package(VTK REQUIRED) INCLUDE_DIRECTORIES(${VTK_INCLUDE_DIRS}) FIND_PACKAGE(GDCM REQUIRED) -include(${GDCM_USE_FILE}) INCLUDE_DIRECTORIES(${GDCM_INCLUDE_DIRS}) SET(SRCS @@ -34,7 +33,7 @@ SET(LINK_LIBS ${LIB_NAME} ${VTK_LIBRARIES} ${ITK_LIBRARIES} - ${GDCM_LIBRARIES} + gdcmMSFF ) TARGET_LINK_LIBRARIES( @@ -45,3 +44,7 @@ TARGET_LINK_LIBRARIES( target_include_directories (${LIB_NAME} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ) + + +set(EXTRA_LIBS ${EXTRA_LIBS} itkDTRrecon) +set(EXTRA_LIBS ${EXTRA_LIBS} PARENT_SCOPE) diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index 59ac4d9..5fb3ad1 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -22,6 +22,7 @@ gfattori 08.11.2021 #include "itkMatrix.h" #include "itkGDCMSeriesFileNames.h" #include "itkOrientImageFilter.h" +#include #include @@ -284,12 +285,12 @@ void itkImageProcessor::SetCustom_ImportTransform(double dTx, } void itkImageProcessor::SetCustom_ProjectionCenterOffsetFixedAxes_IEC( - double dX1, - double dY1, - double dZ1, - double dX2, - double dY2, - double dZ2) + double dX1, + double dY1, + double dZ1, + double dX2, + double dY2, + double dZ2) { typedef itk::Point PointType; @@ -369,8 +370,108 @@ void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2) //TODO UPDATE TO FOLLOW } +int itkImageProcessor::load3DSerieFromFiles( std::vector v_fnames){ -int itkImageProcessor::load3DSerie(const char * pcDirName) + tPatOrientation m_PatOrientation + = tPatOrientation::NotDefined; + + /* Since are not sure DB indexer to sort images, + * we run the IPP sorter here. */ + using IppSorterType = gdcm::IPPSorter; + IppSorterType ZSorter; + + ZSorter.SetComputeZSpacing( true ); + ZSorter.SetZSpacingTolerance( 1e-3 ); + bool b = ZSorter.Sort(v_fnames); + if( !b ) + { + std::cerr << "itkImageProcessor::load3DSerieFromFiles :: Failed to sort files based on IPP." << std::endl; + return EXIT_FAILURE; + } + + /* Some debugging output */ + //std::cout << "Sorting succeeded:" << std::endl; + //ZSorter.Print( std::cout ); + //std::cout << "Found z-spacing:" << std::endl; + //std::cout << ZSorter.GetZSpacing() << std::endl; + + const std::vector & v_sortedFnames = + ZSorter.GetFilenames(); + + try + { + + using ImageIOType = itk::GDCMImageIO; + ImageIOType::Pointer dicomIO = ImageIOType::New(); + + ImageSeriesReaderType3D::Pointer + imageSeriesReader3D = ImageSeriesReaderType3D::New(); + + imageSeriesReader3D->SetImageIO(dicomIO); + imageSeriesReader3D->SetFileNames(v_sortedFnames); + //imageSeriesReader3D->SetNumberOfWorkUnits(8); + imageSeriesReader3D->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt + //imageSeriesReader3D->SetSpacingWarningRelThreshold(); + + imageSeriesReader3D->Update(); + + if(v_fnames.size()>0){ + + gdcm::Reader R; + R.SetFileName(v_sortedFnames.at(0).c_str()); + if(!R.Read()) + { + cerr<<"ERROR: cannot read file: "<SetInput(imageSeriesReader3D->GetOutput()); + caster3D->Update(); + + m_VolumeSourceDupli->SetInputImage(caster3D->GetOutput()); + m_VolumeSourceDupli->Update(); + + caster3D = NULL; + + + } + catch (itk::ExceptionObject & ex) + { + std::cout << ex << std::endl; + return EXIT_FAILURE; + } + + InternalImageType::Pointer m_InputImage = + m_VolumeSourceDupli->GetOutput(); + return + this->fill3DVolumeMeta(m_InputImage, m_PatOrientation); + +} + +int itkImageProcessor::load3DSerieFromFolder(const char * pcDirName) { using NamesGeneratorType = itk::GDCMSeriesFileNames; @@ -402,7 +503,7 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) else { std::cout << "No DICOMs in: " << pcDirName << std::endl; - return EXIT_SUCCESS; + return EXIT_FAILURE; } while (seriesItr != seriesEnd) @@ -467,7 +568,7 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) m_PatOrientation = eImageOrientationType::NotDefined; } - free (sTmpString); + delete [] (sTmpString); } else { return -1; } @@ -493,9 +594,17 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) return EXIT_FAILURE; } - InternalImageType::Pointer m_InputImage = m_VolumeSourceDupli->GetOutput(); + return + this->fill3DVolumeMeta(m_InputImage, m_PatOrientation); + +} + +/** Fill Meta after 3D volume load */ +int itkImageProcessor::fill3DVolumeMeta( + InternalImageType::Pointer m_InputImage, + tPatOrientation m_PatOrientation){ if(m_CTMetaInfo != NULL){ @@ -597,9 +706,9 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) image3DIn = m_3DInputChangeInformationToZero->GetOutput(); - return EXIT_SUCCESS; + } const std::vector itkImageProcessor::GetRTImportOffset() @@ -1176,7 +1285,7 @@ void itkImageProcessor::InitializeProjector() //AAAAAA TODOOOOO CERCA ImageType3D::PointType pFakeIsoLPS = - m_DRTImage1MetaInfo->GetProjectionOriginLPS(); + m_DRTImage1MetaInfo->GetProjectionOriginLPS(); //ImageType3D::PointType pFakeIsoLPS = // m_CTMetaInfo->GetProjectionOriginLPS( @@ -1249,7 +1358,7 @@ void itkImageProcessor::InitializeProjector() if(m_RTMetaInfo == NULL) { ImageType3D::PointType pFakeIsoLPS = - m_DRTImage2MetaInfo->GetProjectionOriginLPS(); + m_DRTImage2MetaInfo->GetProjectionOriginLPS(); //ImageType3D::PointType pFakeIsoLPS = // m_CTMetaInfo->GetProjectionOriginLPS( @@ -1317,17 +1426,17 @@ void itkImageProcessor::InitializeProjector() -// // 2D Image 2 -// interpolator2->SetProjectionAngle( -// dtr * -// m_DRTImage2MetaInfo->GetProjectionAngleLPS() ); -// interpolator2->SetFocalPointToIsocenterDistance( -// m_DRTImage2MetaInfo->GetSCD()); -// interpolator2->SetThreshold( -// m_DRTGeometryMetaInfo->GetIntensityThreshold() -// ); -// interpolator2->SetTransform(transform); -// 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(); @@ -1518,10 +1627,10 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ ImageType3D::PointType NominalIsocenterZero_wZcorrectionLPS; - NominalIsocenterZero_wZcorrectionLPS = - m_CTMetaInfo->GetProjectionOriginLPSZero( - m_DRTGeometryMetaInfo->GetProjectionCenter() - ); + NominalIsocenterZero_wZcorrectionLPS = + m_CTMetaInfo->GetProjectionOriginLPSZero( + m_DRTGeometryMetaInfo->GetProjectionCenter() + ); std::cout<< "///////////////// PGEOM META BEG ///////////////" <GetProjectionOriginLPS() ); - std::cout<<"Bounds1: " + std::cout<<"Bounds1: " <SetSizeWithBounds( vBounds.data(), @@ -1593,8 +1702,8 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ m_DRTImage1MetaInfo->GetCOV(), NominalIsocenter_wZcorrectionLPS,//m_DRTImage1MetaInfo->GetProjectionOriginLPS(), this->CalcProjectionAngleLPS( - m_CTMetaInfo->GetPatientOrientation(), - m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() ), + m_CTMetaInfo->GetPatientOrientation(), + m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() ), //m_DRTGeometryMetaInfo->GetProjectionAngle1IEC(),// m_DRTImage1MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS ) @@ -1605,9 +1714,9 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ m_DRTImage1MetaInfo->SetImageDirectionsLPS( CalcDRTImageDirections( this->CalcProjectionAngleLPS( - m_CTMetaInfo->GetPatientOrientation(), - m_DRTGeometryMetaInfo->GetProjectionAngle1IEC()), - //m_DRTGeometryMetaInfo->GetProjectionAngle1IEC(),//m_DRTImage1MetaInfo->GetProjectionAngleLPS(), + m_CTMetaInfo->GetPatientOrientation(), + m_DRTGeometryMetaInfo->GetProjectionAngle1IEC()), + //m_DRTGeometryMetaInfo->GetProjectionAngle1IEC(),//m_DRTImage1MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS) ); @@ -1648,15 +1757,15 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ m_DRTImage2MetaInfo->SetProjectionOriginLPS( CalibratedIsocenterLPS); -// m_CTMetaInfo->GetProjectionOriginLPS( -// m_DRTGeometryMetaInfo->GetProjectionCenter2()) -// ); + // m_CTMetaInfo->GetProjectionOriginLPS( + // m_DRTGeometryMetaInfo->GetProjectionCenter2()) + // ); m_DRTImage2MetaInfo->SetProjectionOriginLPSZero( CalibratedIsocenterZeroLPS); -// m_CTMetaInfo->GetProjectionOriginLPSZero( -// m_DRTGeometryMetaInfo->GetProjectionCenter2()) -// ); + // m_CTMetaInfo->GetProjectionOriginLPSZero( + // m_DRTGeometryMetaInfo->GetProjectionCenter2()) + // ); vBounds.clear(); vBounds = m_CTMetaInfo->CalculateRegions( @@ -1664,12 +1773,12 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ ); std::cout<<"Bounds2: " - <SetSizeWithBounds( vBounds.data(), @@ -1683,9 +1792,9 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ m_DRTImage2MetaInfo->GetCOV(), NominalIsocenter_wZcorrectionLPS,//m_DRTImage2MetaInfo->GetProjectionOriginLPS(), this->CalcProjectionAngleLPS( - m_CTMetaInfo->GetPatientOrientation(), - m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , - //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), + m_CTMetaInfo->GetPatientOrientation(), + m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , + //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS ) ); @@ -1693,8 +1802,8 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){ m_DRTImage2MetaInfo->SetImageDirectionsLPS( CalcDRTImageDirections( this->CalcProjectionAngleLPS( - m_CTMetaInfo->GetPatientOrientation(), - m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , + m_CTMetaInfo->GetPatientOrientation(), + m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) , //m_DRTGeometryMetaInfo->GetProjectionAngle2IEC(),//m_DRTImage2MetaInfo->GetProjectionAngleLPS(), Std_DRT2LPS) ); @@ -1773,8 +1882,8 @@ void itkImageProcessor::GetProjectionImages(){ ImageType3D::PointType pFakeIsoLPS = m_DRTImage1MetaInfo->GetProjectionOriginLPS(); - //m_CTMetaInfo->GetProjectionOriginLPS( - // m_DRTGeometryMetaInfo->GetProjectionCenter()); + //m_CTMetaInfo->GetProjectionOriginLPS( + // m_DRTGeometryMetaInfo->GetProjectionCenter()); std::cout<<"pFakeIsoLPS: "<GetAngleZ() ); - std::cout<< "itkImageProcessor::GetProjectionImages" <SetCenter( m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()); @@ -1851,65 +1960,65 @@ void itkImageProcessor::GetProjectionImages(){ //CurrTransform->SetComputeZYX(true); //CurrTransform->SetIdentity(); - if(m_RTMetaInfo == NULL) - { + if(m_RTMetaInfo == NULL) + { - ImageType3D::PointType pFakeIsoLPS = - m_DRTImage2MetaInfo->GetProjectionOriginLPS(); - //m_CTMetaInfo->GetProjectionOriginLPS( - // m_DRTGeometryMetaInfo->GetProjectionCenter()); + 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()); - transform2->SetRotation( - CurrTransform->GetAngleX(), - CurrTransform->GetAngleY(), - CurrTransform->GetAngleZ() + CurrTransform = + CalculateInternalTransform( + ZeroPoint, + ZeroPoint, + m_TransformMetaInfo->GetTranslations(), + m_TransformMetaInfo->GetRotations(), + m_DRTImage2MetaInfo->GetProjectionOriginLPS(), + pFakeIsoLPS, + ZeroPoint, + IECtoLPS_Directions ); + std::cout<<"pFakeIsoLPS: "<SetCenter( - m_DRTImage2MetaInfo->GetProjectionOriginLPSZero()); + } 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 ->Print(std::cout); + + 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); @@ -1928,10 +2037,10 @@ void itkImageProcessor::GetProjectionImages(){ imageDRT2In = filter2->GetOutput(); std::cout<< "itkImageProcessor::GetProjectionImages END" <GetInverseTransform()->GetMatrix()<GetInverseTransform()->GetTranslation()<GetInverseTransform()->GetMatrix()<GetInverseTransform()->GetTranslation()< ); + int load2D(const char *); void loadRTPlanInfo(double, double, double, double, double ,double); @@ -151,13 +153,13 @@ public: protected: /** Various pixel types */ - using InternalPixelType = float; - using PixelType2D = short; - using PixelType3D = short; - using OutputPixelType = unsigned char; + using InternalPixelType = float; + using PixelType2D = short; + using PixelType3D = short; + using OutputPixelType = unsigned char; - using ImageType3D = itk::Image; - using InternalImageType = itk::Image; + using ImageType3D = itk::Image; + using InternalImageType = itk::Image; itkImageProcessor(); virtual ~itkImageProcessor(); @@ -167,6 +169,9 @@ private: itkImageProcessor(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented + /** Fill Meta after 3D volume load */ + int fill3DVolumeMeta(InternalImageType::Pointer, + tPatOrientation); /** Image types */ using ImageType2D = itk::Image; diff --git a/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.cxx b/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.cxx index 8853be8..c3bd7bb 100644 --- a/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.cxx +++ b/reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.cxx @@ -148,6 +148,8 @@ int vtkContourTopogramProjectionFilter::RequestData( pp2[1] -= dProjectionLPSOffset [1]; pp2[2] -= dProjectionLPSOffset [2]; +// std::cout << ii << " " << pp2[0] << " " << pp2[1] << " " << pp2[2] << std::endl; + PrjPoints->InsertPoint(ii,pp2); }