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; + };