From d9e7d4fe6d61ff13b50c0f14248737b53a70d5ff Mon Sep 17 00:00:00 2001 From: fattori_g Date: Mon, 7 Feb 2022 18:00:22 +0100 Subject: [PATCH] Rework of projection reference frame. Now we preserve patient LPS until the end and DRT can be directly overlay with Topogram with no further manipulation. Projeciton is in Patient Orientation. Projection angle is set by user in IEC and corrected for patient orientation (HFS,FFS, ..) before internal use Detection of Topogram projection (LAT or PA) is done on load, accounting for patient orientation and image direction cosines --- .../itkDTRrecon/itkImageProcessor.cpp | 972 ++++++++++-------- .../itkDTRrecon/itkImageProcessor.h | 68 +- 2 files changed, 570 insertions(+), 470 deletions(-) diff --git a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp index ee03f7c..13ca1f2 100644 --- a/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp +++ b/reg23Topograms/itkDTRrecon/itkImageProcessor.cpp @@ -21,22 +21,32 @@ gfattori 08.11.2021 #include "itkPermuteAxesImageFilter.h" #include "itkMatrix.h" #include "itkGDCMSeriesFileNames.h" +#include "itkOrientImageFilter.h" #include namespace itk { +/* this is in the end a IEC to HFS... + * but we keep this for ourselves... + */ +static double Standard_DRT2LPS [9] = { + 1,0,0, + 0,0,-1, + 0,1,0 }; -static double coronalElements[9] = { - 1, 0, 0, - 0, 0, 1, - 0,-1, 0 }; -static double sagittalElements[9] = { - 0, 0,-1, - 1, 0, 0, - 0,-1, 0}; + +static double PAElementsIEC[9] = { + 1, 0, 0, + 0, -1, 0, + 0, 0, -1 }; + +static double LATElementsIEC[9] = { + 0, 0, -1, + 0, -1, 0, + -1, 0, 0}; static double HFS2IEC[9] = { 1, 0, 0, @@ -91,7 +101,7 @@ itkImageProcessor::itkImageProcessor() registration->SetInterpolator1(interpolator1); registration->SetInterpolator2(interpolator2); - imageReader2D = ImageReaderType2D::New(); + //imageReader2D = ImageReaderType2D::New(); imageReader3D = ImageReaderType3D::New(); imageSeriesReader3D = ImageSeriesReaderType3D::New(); @@ -129,32 +139,42 @@ itkImageProcessor::itkImageProcessor() using FlipAxesArrayType = FlipFilterType::FlipAxesArrayType; FlipAxesArrayType flipArray; - flipArray[0] = true; - flipArray[1] = true; - flipArray[2] = false; - flipFilter1->SetFlipAxes(flipArray); - //flipFilterLocalizer1->SetFlipAxes(flipArray); flipArray[0] = false; - flipArray[1] = true; + flipArray[1] = false; flipArray[2] = false; +// flipArray[0] = true; +// flipArray[1] = true; +// flipArray[2] = false; + flipFilter1->SetFlipAxes(flipArray); +// flipArray[0] = false; +// flipArray[1] = true; +// flipArray[2] = false; flipFilter2->SetFlipAxes(flipArray); - //flipFilterLocalizer2->SetFlipAxes(flipArray); caster3D = CastFilterType3D::New(); caster2D1 = CastFilterType3D::New(); caster2D2 = CastFilterType3D::New(); - caster2DInput = CastFilterType3D::New(); toVTK2D1 = ITKtoVTKFilterType::New(); toVTK2D2 = ITKtoVTKFilterType::New(); toVTKLocalizer1 = ITKtoVTKFilterType::New(); toVTKLocalizer2 = ITKtoVTKFilterType::New(); - image2D1Rescaler = Input2DRescaleFilterType::New(); - image2D2Rescaler = Input2DRescaleFilterType::New(); - image3DIECFilt = ChangeInformationFilterInputType::New(); + m_3DPatOrientation = eImageOrientationType::NotDefined; + + m_LATSourceDupli = DuplicatorType::New(); + m_PASourceDupli = DuplicatorType::New(); + + + for(int iIdx = 0 ; iIdx < 3; iIdx++){ + for(int jIdx = 0 ; jIdx < 3; jIdx++){ + Std_DRT2LPS.GetVnlMatrix()[iIdx][jIdx] = Standard_DRT2LPS [jIdx+iIdx*3]; + } + } + + } @@ -280,8 +300,6 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) imageSeriesReader3D->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt imageSeriesReader3D->Update(); - image3DIECFilt->SetInput(imageSeriesReader3D->GetOutput()); - //image3DIn = imageSeriesReader3D->GetOutput(); eImageOrientationType currImgOrient; @@ -305,9 +323,12 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::HFS])){ std::cout<<"Image Orientation: HFS"< toIECMatrix; - toIECMatrix.Fill(0.0); - InternalImageType::DirectionType toIECDirection; for(int iIdx = 0 ; iIdx < 3; iIdx++){ for(int jIdx = 0 ; jIdx < 3; jIdx++){ switch (currImgOrient) { case eImageOrientationType::HFS: - toIECDirection.GetVnlMatrix()[iIdx][jIdx] = HFS2IEC[jIdx+iIdx*3]; + LPStoIEC_Directions.GetVnlMatrix()[iIdx][jIdx] = HFS2IEC[jIdx+iIdx*3]; break; case eImageOrientationType::FFS: - toIECDirection.GetVnlMatrix()[iIdx][jIdx] = FFS2IEC[jIdx+iIdx*3]; + LPStoIEC_Directions.GetVnlMatrix()[iIdx][jIdx] = FFS2IEC[jIdx+iIdx*3]; break; default: break; @@ -338,20 +356,12 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) } -// std::cout<< "Orignal Direction: " -// << imageSeriesReader3D->GetOutput()->GetDirection()<GetOutput()->GetDirection() * toIECDirection<SetOutputDirection( -// imageSeriesReader3D->GetOutput()->GetDirection() * toIECDirection -// ); -// image3DIECFilt->ChangeDirectionOn(); + /* Calculate projection angle IEC to LPS */ + projAngle1 = CalcProjectionAngleLPS(currImgOrient, projAngle1_IEC); + projAngle2 = CalcProjectionAngleLPS(currImgOrient, projAngle2_IEC); - image3DIECFilt->Update(); - image3DIn = image3DIECFilt->GetOutput(); + image3DIn = imageSeriesReader3D->GetOutput(); /* before changing origin for projection, let's save the center of volume */ using ImageRegionType3D = ImageType3D::RegionType; using SizeType3D = ImageRegionType3D::SizeType; @@ -363,6 +373,9 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) const itk::Vector resolution3D = image3DIn->GetSpacing(); ImageDirectionType3D imagDirection = image3DIn->GetDirection(); + + + /* calculate image size in 3D Space by fiding the last voxel position */ using VectorType = itk::Vector; VectorType Dest; @@ -373,9 +386,9 @@ int itkImageProcessor::load3DSerie(const char * pcDirName) ImageType3D::PointType LastVoxelPosition = origin3D + (imagDirection * Dest); - image3DCOV[0] = (origin3D[0]+LastVoxelPosition[0])/2; - image3DCOV[1] = (origin3D[1]+LastVoxelPosition[1])/2; - image3DCOV[2] = (origin3D[2]+LastVoxelPosition[2])/2; + image3DCOV[0] = ((origin3D)[0]+LastVoxelPosition[0])/2; + image3DCOV[1] = ((origin3D)[1]+LastVoxelPosition[1])/2; + image3DCOV[2] = ((origin3D)[2]+LastVoxelPosition[2])/2; if(true){ std::cout<<" -------- 3D IMAGE --------"<SetFileName(pcFName); - imageReader3D->Update(); - - image3DIn = imageReader3D->GetOutput(); - - - /* before changing origin for projection, let's save the center of volume */ - using ImageRegionType3D = ImageType3D::RegionType; - using SizeType3D = ImageRegionType3D::SizeType; - using ImageDirectionType3D = ImageType3D::DirectionType; - - ImageRegionType3D region3D = image3DIn->GetBufferedRegion(); - SizeType3D size3D = region3D.GetSize(); - ImageType3D::PointType origin3D = image3DIn->GetOrigin(); - const itk::Vector resolution3D = image3DIn->GetSpacing(); - ImageDirectionType3D imagDirection = image3DIn->GetDirection(); - - /* calculate image size in 3D Space by fiding the last voxel position */ - using VectorType = itk::Vector; - VectorType Dest; - Dest[0]=(size3D[0]-1) * resolution3D[0]; - Dest[1]=(size3D[1]-1) * resolution3D[1]; - Dest[2]=(size3D[2]-1) * resolution3D[2]; - - ImageType3D::PointType LastVoxelPosition = - origin3D + (imagDirection * Dest); - - image3DCOV[0] = (origin3D[0]+LastVoxelPosition[0])/2; - image3DCOV[1] = (origin3D[1]+LastVoxelPosition[1])/2; - image3DCOV[2] = (origin3D[2]+LastVoxelPosition[2])/2; - - if(true){ - std::cout<<" -------- 3D IMAGE --------"<SetOrigin(image3DOrigin); + InternalImageType::DirectionType + DRTDirsIEC; + + DRTDirsIEC.SetIdentity(); + DRTDirsIEC. GetVnlMatrix()[0][0] = cos (pAngleIEC * itk::Math::pi/180); + DRTDirsIEC. GetVnlMatrix()[0][2] = sin (pAngleIEC * itk::Math::pi/180); + DRTDirsIEC. GetVnlMatrix()[2][0] = -sin (pAngleIEC * itk::Math::pi/180); + DRTDirsIEC. GetVnlMatrix()[2][2] = cos (pAngleIEC * itk::Math::pi/180); - caster3D->SetInput(image3DIn); - caster3D->Update(); +// std::cout<<"DRTDirsIEC: "<< DRTDirsIEC < tokens; - /* Intensity window center */ strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0028,0x1050), ds)); for (auto i = strtok(&sTmpString[0], "\\"); i != NULL; i = strtok(NULL, "\\")) @@ -555,13 +584,17 @@ int itkImageProcessor::load2D(const char * pcFName){ std::cout<<"Intensity window center and width: "<< dIntensityWindow[0] <<" "<< dIntensityWindow[1] <SetFileName(pcFName); imageReader2D->SetImageIO(dicomIO); @@ -582,206 +617,163 @@ int itkImageProcessor::load2D(const char * pcFName){ return -1; } - using ImageRegionType3D = ImageType3D::RegionType; - using SizeType3D = ImageRegionType3D::SizeType; + /* identifing a projection to be lateral or posterior requires the interpretation of direction + * and patient orientation */ + using ImageDirectionType3D = ImageType3D::DirectionType; + ImageDirectionType3D LocalizerImagDirectionDCM = + imageReader2D->GetOutput()->GetDirection(); - ImageRegionType3D region3D = imageReader2D->GetOutput()->GetBufferedRegion(); - SizeType3D size3D = region3D.GetSize(); - ImageType3D::PointType origin3D = imageReader2D->GetOutput()->GetOrigin(); - const itk::Vector resolution3D = imageReader2D->GetOutput()->GetSpacing(); - ImageDirectionType3D imagDirection = imageReader2D->GetOutput()->GetDirection(); + std::cout<<"LocalizerImagDirectionDCM " <; - VectorType Dest; - Dest[0]=(size3D[0]-1) * resolution3D[0]; - Dest[1]=(size3D[1]-1) * resolution3D[1]; - Dest[2]=(size3D[2]-1) * resolution3D[2]; - ImageType3D::PointType LastVoxelPosition = - origin3D + (imagDirection * Dest); + /* calculate image orientation with respect to fixed IEC */ + InternalImageType::DirectionType LocalizerImagDirectionIEC = + toIECDirection * LocalizerImagDirectionDCM; - if(false){ - std::cout<<" -------- LOAD OF 2D IMAGE --------"<GetOutput(); - - // for(int iIdx =0 ; iIdx <100; iIdx++){ - // image_tmp1.GetPointer()->GetBufferPointer()[iIdx] = 2000; - // } - - /* cast type to internal processing type */ - caster2DInput->SetInput(image_tmp1); - caster2DInput->Update(); + double dSumLAT = 0; + dSumLAT= + (LocalizerImagDirectionIEC[0][2] * LATElementsIEC[2]) + + (LocalizerImagDirectionIEC[1][2] * LATElementsIEC[5]) + + (LocalizerImagDirectionIEC[2][2] * LATElementsIEC[8]); - using DuplicatorType = itk::ImageDuplicator; - DuplicatorType::Pointer duplicator = DuplicatorType::New(); - duplicator->SetInputImage(caster2DInput->GetOutput()); - duplicator->Update(); + tProjOrientationType + currProjOrientation = NA; - ImageType3D::PointType image3DOrigin; - image3DOrigin[0]=image3DCOV[0]; - image3DOrigin[1]=image3DCOV[1]; - image3DOrigin[2]=image3DCOV[2]; + /* dot shoudl be 1 or -1 if the normal is parallel */ + if(dSumPA == 1 || dSumPA == -1){ + currProjOrientation = PA; + std::cout<<"It's a PA image"<SetInput(imageReader2D->GetOutput()); + caster2DInput->Update(); - image2D1In = duplicator->GetOutput(); + InternalImageType::Pointer MyLocalizerImage = + caster2DInput->GetOutput(); - const InternalImageType::DirectionType image2D1InDirection = image2D1In->GetDirection(); - // std::cout<<"image2D1InDirection "<GetDirection().GetTranspose(); - This is the distance between the topogram origin and the 3DVolume center (which - are given in DICOM FOR) mapped into the topogram frame of reference + image2Dcenter_3DCOV = - + ( + image2D1InDirectionInv * (image3DOrigin - MyLocalizerImage->GetOrigin()) + ); - */ - ImageType3D::PointType image2Dcenter_3DCOV; - InternalImageType::DirectionType image2D1InDirectionInv; - image2D1InDirectionInv.GetVnlMatrix() = image2D1InDirection.GetTranspose(); + /* + We have to set the z coordinate to the focal distance, it would be 0 otherwise... + */ + image2Dcenter_3DCOV[2] = -scd; - image2Dcenter_3DCOV = - - ( - image2D1InDirectionInv * (image3DOrigin - image2D1In->GetOrigin()) - ); + ChangeInformationFilterType::Pointer filter = ChangeInformationFilterType::New(); + filter->SetInput(MyLocalizerImage); - /* - We have to set the z coordinate to the focal distance, it would be 0 otherwise... - */ - image2Dcenter_3DCOV[2] = -scd; + const InternalImageType::DirectionType newDirection = + MyLocalizerImage->GetDirection() * image2D1InDirectionInv; +// const InternalImageType::DirectionType newDirection = +// toIECDirection * MyLocalizerImage->GetDirection() ; - ChangeInformationFilterType::Pointer filter = ChangeInformationFilterType::New(); - filter->SetInput(image2D1In); + //filter->SetOutputDirection(newDirection); + //filter->ChangeDirectionOn(); + //filter->SetOutputOrigin(image2Dcenter_3DCOV); + //filter->ChangeOriginOn(); + try + { + filter->UpdateOutputInformation(); + } + catch (itk::ExceptionObject & error) + { + std::cerr << "Error: " << error << std::endl; + return EXIT_FAILURE; + } - const InternalImageType::DirectionType newDirection = image2D1InDirection * image2D1InDirectionInv;//direction * rotation.GetMatrix(); - - filter->SetOutputDirection(newDirection); - filter->ChangeDirectionOn(); - - filter->SetOutputOrigin(image2Dcenter_3DCOV); - filter->ChangeOriginOn(); - - try - { - filter->UpdateOutputInformation(); - } - catch (itk::ExceptionObject & error) - { - std::cerr << "Error: " << error << std::endl; - return EXIT_FAILURE; - } + filter->Update(); - memcpy(image1IntensityWindow,dIntensityWindow, 2*sizeof(double)); + InternalImageType::Pointer m_Image; + double* m_ImageIntensity; + FlipFilterType::Pointer m_FlipFilter; + DuplicatorType::Pointer m_Duplicator; + bool* m_ImLoaded; + + switch (currProjOrientation) { + case eProjectionOrientationType::PA: + + m_Image = image2D1In.GetPointer(); + m_ImageIntensity = image1IntensityWindow; + m_FlipFilter = flipFilterLocalizer1.GetPointer(); + m_Duplicator = m_PASourceDupli.GetPointer(); + m_ImLoaded = &image2D1Loaded; + + break; + + case eProjectionOrientationType::LAT: + + m_Image = image2D2In.GetPointer(); + m_ImageIntensity = image2IntensityWindow; + m_FlipFilter = flipFilterLocalizer2.GetPointer(); + m_Duplicator = m_LATSourceDupli.GetPointer(); + m_ImLoaded = &image2D2Loaded; + + break; + + default: + return -1; + break; + } - flipFilterLocalizer1->SetInput(filter->GetOutput()); - flipFilterLocalizer1->Update(); + memcpy(m_ImageIntensity,dIntensityWindow, 2*sizeof(double)); - image2D1Rescaler->SetInput(flipFilterLocalizer1->GetOutput()); - image2D1Rescaler->SetOutputMinimum(0); - image2D1Rescaler->SetOutputMaximum(255); + m_Duplicator->SetInputImage(filter->GetOutput()); + m_Duplicator->Update(); + m_Image = m_Duplicator->GetOutput(); + + m_FlipFilter->SetInput(m_Image); - image2D1Loaded = true; - return 1; - - } else if (dSumSagittal < 1e-10){ - - //std::cout<< "Sagittal" <GetOutput(); - - const InternalImageType::DirectionType image2D2InDirection = image2D2In->GetDirection(); - - /* center topogram at the COV of 3D image - this works only if user does not set custom isocenter - to be fixed. - - This is the distance between the topogram origin and the 3DVolume center (which - are given in DICOM FOR) mapped into the topogram frame of reference - - */ - ImageType3D::PointType image2Dcenter_3DCOV; - InternalImageType::DirectionType image2D2InDirectionInv; - image2D2InDirectionInv.GetVnlMatrix() = image2D2InDirection.GetTranspose(); - - image2Dcenter_3DCOV = - - ( - image2D2InDirectionInv * (image3DOrigin - image2D2In->GetOrigin()) - ); - - /* - We have to set the z coordinate to the focal distance, it would be 0 otherwise... - */ - - image2Dcenter_3DCOV[2] = -scd; - - - using FilterType = itk::ChangeInformationImageFilter; - FilterType::Pointer filter = FilterType::New(); - filter->SetInput(image2D2In); - - const InternalImageType::DirectionType newDirection = image2D2InDirection * image2D2InDirectionInv;//image2D2InDirectionInv; - filter->SetOutputDirection(newDirection); - filter->ChangeDirectionOn(); - - filter->SetOutputOrigin(image2Dcenter_3DCOV); - filter->ChangeOriginOn(); - - try - { - filter->UpdateOutputInformation(); - } - catch (itk::ExceptionObject & error) - { - std::cerr << "Error: " << error << std::endl; - return EXIT_FAILURE; - } - - memcpy(image2IntensityWindow,dIntensityWindow, 2*sizeof(double)); - - - flipFilterLocalizer2->SetInput(filter->GetOutput()); - flipFilterLocalizer2->Update(); - - image2D2Rescaler->SetInput(flipFilterLocalizer2->GetOutput()); - image2D2Rescaler->SetOutputMinimum(0); - image2D2Rescaler->SetOutputMaximum(255); - - - image2D2Loaded = true; - return 2; - } else { - std::cout<< "Not Accepted" < + +itkImageProcessor::InternalImageType::DirectionType + itkImageProcessor::CalcDRTImageDirections(double angle){ + + + InternalImageType::DirectionType + DRTDirs; + + DRTDirs.SetIdentity(); + DRTDirs. GetVnlMatrix()[0][0] = cos (angle*itk::Math::pi/180); + DRTDirs. GetVnlMatrix()[0][1] = -sin (angle*itk::Math::pi/180); + DRTDirs. GetVnlMatrix()[1][0] = sin (angle*itk::Math::pi/180); + DRTDirs. GetVnlMatrix()[1][1] = cos (angle*itk::Math::pi/180); + + DRTDirs *= Std_DRT2LPS; + + return + DRTDirs; + +} + + +/* + Nobody should go through all this... + To be very careful editing ... + +*/ +itkImageProcessor::ImageType3D::PointType + itkImageProcessor::CalcDRTImageOrigin(itkImageProcessor::InternalImageType::Pointer m_Image, + itkImageProcessor::ImageType3D::PointType m_VolCOOV, + double dAngle + ){ + + using ImageRegionType3D = itkImageProcessor::ImageType3D::RegionType; + using SizeType3D = ImageRegionType3D::SizeType; + using ImageDirectionType3D = itkImageProcessor::ImageType3D::DirectionType; + using VectorType = itk::Vector; + + ImageRegionType3D region3D; + SizeType3D size3D; + itk::Vector resolution3D; + + + region3D = m_Image ->GetBufferedRegion(); + size3D = region3D.GetSize(); + resolution3D = m_Image ->GetSpacing(); + + VectorType Dest; + Dest[0]=(size3D[0]-1) * resolution3D[0]; + Dest[1]=(size3D[1]-1) * resolution3D[1]; + Dest[2]=(size3D[2]-1) * resolution3D[2]; + + itkImageProcessor::ImageType3D::PointType LastVoxelPosition = + m_Image ->GetOrigin() + + (m_Image ->GetDirection() * Dest); + + itkImageProcessor::ImageType3D::PointType DRTCOV; + DRTCOV[0] = (( m_Image ->GetOrigin())[0]+LastVoxelPosition[0])/2; + DRTCOV[1] = (( m_Image ->GetOrigin())[1]+LastVoxelPosition[1])/2; + DRTCOV[2] = (( m_Image ->GetOrigin())[2]+LastVoxelPosition[2])/2; + + itkImageProcessor::InternalImageType::DirectionType DRT2LPS + = this->CalcDRTImageDirections(dAngle); + + ImageType3D::PointType NewOrigin = + m_VolCOOV + DRT2LPS * (m_Image ->GetOrigin() - DRTCOV); + + return + NewOrigin; +} + void itkImageProcessor::GetProjectionImages(){ @@ -985,16 +1051,85 @@ void itkImageProcessor::GetProjectionImages(){ } + resampleFilter1->Update(); + resampleFilter2->Update(); - // As explained before, the computed projection is upsided-down. - // Here we use a FilpImageFilter to flip the images in y-direction. - flipFilter1->SetInput(resampleFilter1->GetOutput()); - flipFilter2->SetInput(resampleFilter2->GetOutput()); + /* here probably we should move everything back to dicom */ + + /* Projection images are in a fancy reference system... + * With some knowledge (Projection angle, 3DCOV, etc) they can be + * brought back into Patient reference so that can directly + * overlay with the loaded Localizer images. + */ + + ChangeInformationFilterType::Pointer filter1 = + ChangeInformationFilterType::New(); + + ChangeInformationFilterType::Pointer filter2 = + ChangeInformationFilterType::New(); + + + /* First of all we need the center of the Projection images in its reference frame */ + + ImageType3D::PointType image3DOrigin; + image3DOrigin[0]=image3DCOV[0]; + image3DOrigin[1]=image3DCOV[1]; + image3DOrigin[2]=image3DCOV[2]; + + resampleFilter1->Update(); + + ImageType3D::PointType NewOrigin = + this->CalcDRTImageOrigin(resampleFilter1->GetOutput(), + image3DOrigin, + projAngle1); + InternalImageType::DirectionType DRT2LPS1 = + this->CalcDRTImageDirections(projAngle1); + + filter1->SetInput( resampleFilter1->GetOutput() ); + filter1->SetOutputDirection(DRT2LPS1 );//IECtoLPS_Directions); + filter1->ChangeDirectionOn(); + filter1->SetOutputOrigin( NewOrigin ); + filter1->ChangeOriginOn(); + filter1->UpdateOutputInformation(); + + + + InternalImageType::DirectionType LATtoHFS = + this->CalcDRTImageDirections(projAngle2); + + + /* Again.. */ + resampleFilter2->Update(); + + NewOrigin = + this->CalcDRTImageOrigin(resampleFilter2->GetOutput(), + image3DOrigin, + projAngle2); + + DRT2LPS1 = + this->CalcDRTImageDirections(projAngle2); + + + + filter2->SetInput( resampleFilter2 ->GetOutput() ); + filter2->SetOutputDirection(DRT2LPS1 ); + filter2->ChangeDirectionOn(); + filter2->SetOutputOrigin(NewOrigin ); + filter2->ChangeOriginOn(); + filter2->UpdateOutputInformation(); + + filter1->Update(); + filter2->Update(); + +// As explained before, the computed projection is upsided-down. +// Here we use a FilpImageFilter to flip the images in y-direction. + /* no flip required because now we handle properly the image directions. */ + flipFilter1->SetInput(filter1->GetOutput()); + flipFilter2->SetInput(filter2->GetOutput()); flipFilter1->Update(); flipFilter2->Update(); - // std::cout<< "itkImageProcessor::GetProjectionImages DONE" <; - -//// RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New(); -//// rescaler1->SetOutputMinimum(0); -//// rescaler1->SetOutputMaximum(255); -//// rescaler1->SetInput(flipFilterLocalizer1->GetOutput()); - -//// rescaler1->Update(); -// using ImageRegionType3D = ImageType3D::RegionType; -// using SizeType3D = ImageRegionType3D::SizeType; -// using ImageDirectionType3D = ImageType3D::DirectionType; - -// ImageRegionType3D region3D = image2D1In->GetBufferedRegion(); -// SizeType3D size3D = region3D.GetSize(); -// ImageType3D::PointType origin3D = image2D1In->GetOrigin(); -// const itk::Vector resolution3D = image2D1In->GetSpacing(); -// ImageDirectionType3D imagDirection = image2D1In->GetDirection(); - -// /* calculate image size in 3D Space */ -// using VectorType = itk::Vector; -// VectorType Dest; -// Dest[0]=(size3D[0]-1) * resolution3D[0]; -// Dest[1]=(size3D[1]-1) * resolution3D[1]; -// Dest[2]=(size3D[2]-1) * resolution3D[2]; - -// ImageType3D::PointType LastVoxelPosition = -// origin3D + (imagDirection * Dest); - -// std::cout<<" -------- Localizer 1 ITK --------"<SetInput("image2D1In",image2D1In); -// toVTKLocalizer1->Update(); - -// double* dBounds = toVTKLocalizer1->GetOutput()->GetBounds(); - -// std::cout<< "-------- Localizer 1 VTK --------" <GetOutput(); -//} - -//vtkImageData* itkImageProcessor::GetLocalizer2VTK() -//{ -//// // Rescale the intensity of the projection images to 0-255 for output. -//// using RescaleFilterType = itk::RescaleIntensityImageFilter; - -//// RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New(); -//// rescaler2->SetOutputMinimum(0); -//// rescaler2->SetOutputMaximum(255); -//// rescaler2->SetInput(flipFilterLocalizer2->GetOutput()); - -//// rescaler2->Update(); -// using ImageRegionType3D = ImageType3D::RegionType; -// using SizeType3D = ImageRegionType3D::SizeType; -// using ImageDirectionType3D = ImageType3D::DirectionType; - -// ImageRegionType3D region3D = image2D2In->GetBufferedRegion(); -// SizeType3D size3D = region3D.GetSize(); -// ImageType3D::PointType origin3D = image2D2In->GetOrigin(); -// const itk::Vector resolution3D = image2D2In->GetSpacing(); -// ImageDirectionType3D imagDirection = image2D2In->GetDirection(); - -// /* calculate image size in 3D Space */ -// using VectorType = itk::Vector; -// VectorType Dest; -// Dest[0]=(size3D[0]-1) * resolution3D[0]; -// Dest[1]=(size3D[1]-1) * resolution3D[1]; -// Dest[2]=(size3D[2]-1) * resolution3D[2]; - -// ImageType3D::PointType LastVoxelPosition = -// origin3D + (imagDirection * Dest); - -// std::cout<<" -------- Localizer 2 ITK --------"<SetInput("image2D2In",image2D2In); -// toVTKLocalizer2->Update(); - -// double* dBounds = toVTKLocalizer2->GetOutput()->GetBounds(); - -// std::cout<< "-------- Localizer 2 VTK --------" <GetOutput(); -//} vtkImageData* itkImageProcessor::GetLocalizer1VTK() { @@ -1174,43 +1203,48 @@ vtkImageData* itkImageProcessor::GetLocalizer1VTK() rescaler1->SetInput(flipFilterLocalizer1->GetOutput()); rescaler1->Update(); - using ImageRegionType3D = ImageType3D::RegionType; - using SizeType3D = ImageRegionType3D::SizeType; - using ImageDirectionType3D = ImageType3D::DirectionType; - - ImageRegionType3D region3D = rescaler1->GetOutput()->GetBufferedRegion(); - SizeType3D size3D = region3D.GetSize(); - ImageType3D::PointType origin3D = rescaler1->GetOutput()->GetOrigin(); - const itk::Vector resolution3D = rescaler1->GetOutput()->GetSpacing(); - ImageDirectionType3D imagDirection = rescaler1->GetOutput()->GetDirection(); - - /* calculate image size in 3D Space */ - using VectorType = itk::Vector; - VectorType Dest; - Dest[0]=(size3D[0]-1) * resolution3D[0]; - Dest[1]=(size3D[1]-1) * resolution3D[1]; - Dest[2]=(size3D[2]-1) * resolution3D[2]; - - ImageType3D::PointType LastVoxelPosition = - origin3D + (imagDirection * Dest); - - std::cout<<" -------- Localizer 1 ITK --------"<SetInput(rescaler1->GetOutput()); toVTKLocalizer1->Update(); - double* dBounds = toVTKLocalizer1->GetOutput()->GetBounds(); + if(true) { + using ImageRegionType3D = ImageType3D::RegionType; + using SizeType3D = ImageRegionType3D::SizeType; + using ImageDirectionType3D = ImageType3D::DirectionType; - std::cout<< "-------- Localizer 1 VTK --------" <GetOutput()->GetBufferedRegion(); + SizeType3D size3D = region3D.GetSize(); + ImageType3D::PointType origin3D = rescaler1->GetOutput()->GetOrigin(); + const itk::Vector resolution3D = rescaler1->GetOutput()->GetSpacing(); + ImageDirectionType3D imagDirection = rescaler1->GetOutput()->GetDirection(); + + /* calculate image size in 3D Space */ + using VectorType = itk::Vector; + VectorType Dest; + Dest[0]=(size3D[0]-1) * resolution3D[0]; + Dest[1]=(size3D[1]-1) * resolution3D[1]; + Dest[2]=(size3D[2]-1) * resolution3D[2]; + + ImageType3D::PointType LastVoxelPosition = + origin3D + (imagDirection * Dest); + + std::cout<<" -------- Localizer 1 ITK --------"<GetOutput()->GetBounds(); + + std::cout<< "-------- Localizer 1 VTK --------" <GetOutput(); @@ -1227,45 +1261,51 @@ vtkImageData* itkImageProcessor::GetLocalizer2VTK() rescaler2->SetInput(flipFilterLocalizer2->GetOutput()); rescaler2->Update(); - using ImageRegionType3D = ImageType3D::RegionType; - using SizeType3D = ImageRegionType3D::SizeType; - using ImageDirectionType3D = ImageType3D::DirectionType; - - ImageRegionType3D region3D = rescaler2->GetOutput()->GetBufferedRegion(); - SizeType3D size3D = region3D.GetSize(); - ImageType3D::PointType origin3D = rescaler2->GetOutput()->GetOrigin(); - const itk::Vector resolution3D = rescaler2->GetOutput()->GetSpacing(); - ImageDirectionType3D imagDirection = rescaler2->GetOutput()->GetDirection(); - - /* calculate image size in 3D Space */ - using VectorType = itk::Vector; - VectorType Dest; - Dest[0]=(size3D[0]-1) * resolution3D[0]; - Dest[1]=(size3D[1]-1) * resolution3D[1]; - Dest[2]=(size3D[2]-1) * resolution3D[2]; - - ImageType3D::PointType LastVoxelPosition = - origin3D + (imagDirection * Dest); - - std::cout<<" -------- Localizer 2 ITK --------"<SetInput(rescaler2->GetOutput()); toVTKLocalizer2->Update(); - double* dBounds = toVTKLocalizer2->GetOutput()->GetBounds(); - std::cout<< "-------- Localizer 2 VTK --------" <GetOutput()->GetBufferedRegion(); + SizeType3D size3D = region3D.GetSize(); + ImageType3D::PointType origin3D = rescaler2->GetOutput()->GetOrigin(); + const itk::Vector resolution3D = rescaler2->GetOutput()->GetSpacing(); + ImageDirectionType3D imagDirection = rescaler2->GetOutput()->GetDirection(); + + /* calculate image size in 3D Space */ + using VectorType = itk::Vector; + VectorType Dest; + Dest[0]=(size3D[0]-1) * resolution3D[0]; + Dest[1]=(size3D[1]-1) * resolution3D[1]; + Dest[2]=(size3D[2]-1) * resolution3D[2]; + + ImageType3D::PointType LastVoxelPosition = + origin3D + (imagDirection * Dest); + + std::cout<<" -------- Localizer 2 ITK --------"<GetOutput()->GetBounds(); + std::cout<< "-------- Localizer 2 VTK --------" <GetOutput(); } @@ -1281,6 +1321,11 @@ vtkImageData* itkImageProcessor::GetProjection1VTK() rescaler1->SetInput(flipFilter1->GetOutput()); rescaler1->Update(); + + toVTK2D1->SetInput(rescaler1->GetOutput()); + toVTK2D1->Update(); + + using ImageRegionType3D = ImageType3D::RegionType; using SizeType3D = ImageRegionType3D::SizeType; using ImageDirectionType3D = ImageType3D::DirectionType; @@ -1309,10 +1354,6 @@ vtkImageData* itkImageProcessor::GetProjection1VTK() std::cout<<"Last Voxel position: "<SetInput(rescaler1->GetOutput()); - toVTK2D1->Update(); - double* dBounds = toVTK2D1->GetOutput()->GetBounds(); std::cout<< "-------- Proj 1 --------" <SetInput(flipFilter2->GetOutput()); rescaler2->Update(); + + rescaler2->Print(std::cout); + using ImageRegionType3D = ImageType3D::RegionType; using SizeType3D = ImageRegionType3D::SizeType; using ImageDirectionType3D = ImageType3D::DirectionType; @@ -1369,7 +1413,7 @@ vtkImageData* itkImageProcessor::GetProjection2VTK() double* dBounds = toVTK2D2->GetOutput()->GetBounds(); - std::cout<< "-------- Proj 1 --------" <; + using InternalImageType = itk::Image; + typedef enum eProjectionOrientationType + {NA = 0, LAT=2, PA=1} + tProjOrientationType; + /** Enum type for Orientation */ + typedef enum eImageOrientationType + { NotDefined = 0, HFS = 1, FFS = 2} + tPatOrientation; protected: itkImageProcessor(); virtual ~itkImageProcessor(); void PrintSelf(std::ostream& os, Indent indent) const; + + + private: itkImageProcessor(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented + void Finalize(); -/** Various pixel types */ - using InternalPixelType = float; - using PixelType2D = short; - using PixelType3D = short; - using OutputPixelType = unsigned char; /** Image types */ using ImageType2D = itk::Image; - using ImageType3D = itk::Image; + using OutputImageType = itk::Image; /** The following lines define each of the components used in the registration: The transform, optimizer, metric, interpolator and the registration method itself. */ - using InternalImageType = itk::Image; using TransformType = itk::Euler3DTransform; using OptimizerType = itk::PowellOptimizer; using MetricType = itk::NormalizedCorrelationTwoImageToOneImageMetric; @@ -190,6 +205,7 @@ private: using ImageReaderType3D = itk::ImageFileReader; using ImageSeriesReaderType3D = itk::ImageSeriesReader; using ImageIOType = itk::GDCMImageIO; + using DuplicatorType = itk::ImageDuplicator; /** widely used filters: flip, cast, rescale, ... */ using FlipFilterType = itk::FlipImageFilter; @@ -197,12 +213,12 @@ private: using CastFilterType3D = itk::CastImageFilter; using ChangeInformationFilterType = itk::ChangeInformationImageFilter; using ChangeInformationFilterInputType = itk::ChangeInformationImageFilter; + using ChangeInformationFilterInput2DType = itk::ChangeInformationImageFilter; /** ITK to VTK filter */ using ITKtoVTKFilterType = itk::ImageToVTKImageFilter; - /** Enum type for Orientation */ - enum eImageOrientationType { HFS = 0, FFS = 1}; + MetricType::Pointer metric; @@ -212,7 +228,7 @@ private: InterpolatorType::Pointer interpolator2; RegistrationType::Pointer registration; - ImageReaderType2D::Pointer imageReader2D; +// ImageReaderType2D::Pointer imageReader2D; ImageReaderType3D::Pointer imageReader3D; ImageSeriesReaderType3D::Pointer imageSeriesReader3D; @@ -226,7 +242,7 @@ private: CastFilterType3D::Pointer caster3D; CastFilterType3D::Pointer caster2D1; CastFilterType3D::Pointer caster2D2; - CastFilterType3D::Pointer caster2DInput; + TransformType::InputPointType isocenter3D; @@ -238,15 +254,33 @@ private: ImageType3D::Pointer image3DIn; InternalImageType::Pointer image2D1In; InternalImageType::Pointer image2D2In; + ImageType2D::Pointer image2DInGeneric; - Input2DRescaleFilterType::Pointer image2D1Rescaler; - Input2DRescaleFilterType::Pointer image2D2Rescaler; InternalImageType::Pointer TMPImg1; InternalImageType::Pointer TMPImg2; + DuplicatorType::Pointer m_LATSourceDupli; + DuplicatorType::Pointer m_PASourceDupli; - ChangeInformationFilterInputType::Pointer image3DIECFilt; + tPatOrientation + m_3DPatOrientation; + InternalImageType::DirectionType + LPStoIEC_Directions; + + + InternalImageType::DirectionType + CalcDRTImageDirections(double angle); + + ImageType3D::PointType + CalcDRTImageOrigin(InternalImageType::Pointer m_Image, + ImageType3D::PointType m_VolCOOV, + double dAngle + ); + double + CalcProjectionAngleLPS( + tPatOrientation pOrient, + double pAngleIEC); double image1res[2], image2res[2]; double image1center[2], image2center[2]; @@ -255,6 +289,7 @@ private: double scd; double threshold; double projAngle1, projAngle2; + double projAngle1_IEC, projAngle2_IEC; double TZero[3], RZero[3]; double image3DCOV[3]; @@ -264,6 +299,9 @@ private: double image1IntensityWindow[2], image2IntensityWindow[2]; + + InternalImageType::DirectionType Std_DRT2LPS; + };