fattori_g/reg23Topograms
0
0
Fork 0
Code Issues 1 Pull Requests Actions Packages Projects Releases Wiki Activity

all sparce variables removed.

Browse Source 
Meta for transform implemented.
Runs ok with volume, rt, struct
This commit is contained in:
fattori_g
2022-02-18 00:05:44 +01:00
parent 93b51ee963
commit 34796f0d13
4 changed files with 271 additions and 442 deletions
Download Patch File Download Diff File Expand all files Collapse all files

33
reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp
View File

@ -314,6 +314,39 @@ RTGeometryMetaInformation
TransformMetaInformation ::
TransformMetaInformation (){
m_Translations.Fill(0.);
m_Rotations.Fill(0.);
m_ReferenceTransform.Fill(0.);
m_CurrentTransform.Fill(0.);
}
void
TransformMetaInformation
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
TransformMetaInformation
::~TransformMetaInformation ()
{
}
}

56
reg23Topograms/itkDTRrecon/DRTMetaInformation.h
View File

@ -217,7 +217,6 @@ public:
itkSetMacro(ImageDirections,DirectionType);
itkGetMacro(ImageDirections,DirectionType);
//itkSetMacro(LPS2IECDirections,DirectionType);
itkGetMacro(LPS2IECDirections,DirectionType);
PointType GetCOV();
@ -328,7 +327,10 @@ protected:
m_DRT1Spacing,
m_DRT2Spacing;
PointType
/* Transform between IEC Support and
* IEC Scanner frame of reference */
m_IECS2IECScannerT,
m_IECS2IECScannerR;
@ -353,7 +355,6 @@ private:
class RTGeometryMetaInformation :
public itk::Object{
@ -384,8 +385,12 @@ public:
protected:
PointType
/* RT Plan isocenter in LPS coordinates */
m_IsocenterLPS,
/* RT Plan couch SETUP offset corresponding to
* the RT Plan isocenter */
m_IsocenterIECS;
/** Default Constructor **/
@ -401,6 +406,53 @@ private:
void operator=(const Self&);
};
class TransformMetaInformation :
public itk::Object{
public:
/** standard typedefs **/
typedef TransformMetaInformation Self;
typedef itk::Object Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::Point<double, 3> PointType;
typedef itk::Matrix<double,3,3> TransformMatrixType;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(TransformMetaInformation, itk::Object);
/** object information streaming **/
void PrintSelf(std::ostream& os, itk::Indent indent) const;
protected:
PointType
m_Translations,
m_Rotations;
TransformMatrixType
m_ReferenceTransform,
m_CurrentTransform;
/** Default Constructor **/
TransformMetaInformation ();
/** Default Destructor **/
virtual ~TransformMetaInformation ();
private:
/** purposely not implemented **/
TransformMetaInformation (const Self&);
/** purposely not implemented **/
void operator=(const Self&);
};
}
#endif

561
reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
View File

@ -25,12 +25,6 @@ gfattori 08.11.2021
#include <string>
/* Some parameters that need to be in settings or read from files */
/* Settings */
#define CT_DEFAULT_TABLEHEIGHT 175
namespace itk
@ -44,8 +38,6 @@ static double Standard_DRT2LPS [9] = {
0,0,-1,
0,1,0 };
static double PAElementsIEC[9] = {
1, 0, 0,
0, -1, 0,
@ -97,31 +89,14 @@ itkImageProcessor::itkImageProcessor()
interpolator1 = InterpolatorType::New();
interpolator2 = InterpolatorType::New();
// image1res[0]=image2res[1]= 2.;
// image1Size[0] = image1Size[1] = 512;
// image2Size[0] = image2Size[1] = 512;
TZero[0]=TZero[1]=TZero[2]=0.;
RZero[0]=RZero[1]=RZero[2]=0.;
customized_2DRES = false;
customized_2DSIZE = false;
customized_ImportTransform = false;
customized_ProjectionCenter = false;
customized_RTCouchSetup = false;
customized_RTIsocenter = false;
image2D1Loaded= false;
image2D2Loaded= false;
image3DLoaded= false;
toVTK2D1 = ITKtoVTKFilterType::New();
toVTK2D2 = ITKtoVTKFilterType::New();
toVTKLocalizer1 = ITKtoVTKFilterType::New();
toVTKLocalizer2 = ITKtoVTKFilterType::New();
// m_3DPatOrientation = eImageOrientationType::NotDefined;
m_LATSourceDupli = DuplicatorType::New();
m_PASourceDupli = DuplicatorType::New();
@ -137,12 +112,6 @@ itkImageProcessor::itkImageProcessor()
m_3DInputChangeInformationToZero =
ChangeInformationFilterType::New();
// IEC2DCMMapT.Fill(0.);
// IEC2DCMMapR.Fill(0.);
// rtIsocenterLPS.Fill(0.);
// rtCouchOffset.Fill(0.);
//ProjectionCenterFixedAxes.Fill(0.);
m_Projection1VTKTransform = vtkMatrix4x4::New();
m_Projection1VTKTransform->Identity();
@ -162,6 +131,11 @@ itkImageProcessor::itkImageProcessor()
m_RTMetaInfo = NULL;
m_DRTGeometryMetaInfo = NULL;
m_TransformMetaInfo = NULL;
/* Initialiser the projection geoemtry with defaults */
m_DRTGeometryMetaInfo
= DRTProjectionGeometryImageMetaInformation::New();
@ -195,9 +169,6 @@ itkImageProcessor::itkImageProcessor()
m_DRTGeometryMetaInfo->SetIECS2IECScannerT(Point3D);
m_DRTGeometryMetaInfo->SetIECS2IECScannerR(Point3D);
}
itkImageProcessor::~itkImageProcessor()
@ -213,7 +184,6 @@ void itkImageProcessor::PrintSelf(std::ostream& os, Indent indent) const
void itkImageProcessor::SetIntensityThreshold(double dT)
{
m_DRTGeometryMetaInfo->SetIntensityThreshold(dT);
// d_intensityThreshold = dT;
}
void itkImageProcessor::SetSCD(double dDist)
@ -233,13 +203,6 @@ void itkImageProcessor::SetCustom_ImportTransform(double dTx,
double dRy,
double dRz)
{
// IEC2DCMMapT[0] = dTx;
// IEC2DCMMapT[1] = dTy;
// IEC2DCMMapT[2] = dTz;
// IEC2DCMMapR[0] = dRx;
// IEC2DCMMapR[1] = dRy;
// IEC2DCMMapR[2] = dRz;
customized_ImportTransform = true;
ImageType3D::PointType
Punto;
@ -257,29 +220,9 @@ void itkImageProcessor::SetCustom_ImportTransform(double dTx,
}
//void itkImageProcessor::SetCustom_RTCouchSetup_IEC(double dLAT,
// double dVRT,
// double dLNG)
//{
// rtCouchOffset[0] = dLAT;
// rtCouchOffset[1] = dVRT;
// rtCouchOffset[2] = dLNG;
// customized_RTCouchSetup = true;
//}
//void itkImageProcessor::SetCustom_RTIsocenter_LPS(double dX,
// double dY,
// double dZ)
//{
// rtIsocenterLPS[0] = dX;
// rtIsocenterLPS[1] = dY;
// rtIsocenterLPS[2] = dZ;
// customized_RTIsocenter = true;
//}
void itkImageProcessor::SetCustom_ProjectionCenterFixedAxes_IEC(double dX,
double dY,
double dZ)
double dY,
double dZ)
{
typedef itk::Point<double, 3> PointType;
@ -291,16 +234,11 @@ void itkImageProcessor::SetCustom_ProjectionCenterFixedAxes_IEC(double dX,
m_DRTGeometryMetaInfo->SetProjectionCenter(m_point);
customized_ProjectionCenter = true;
}
void itkImageProcessor::SetCustom_2Dres(double nX1,double nY1,double nX2,double nY2)
{
// image1res[0] = nX1;
// image1res[1] = nY1;
// image2res[0] = nX2;
// image2res[1] = nY2;
customized_2DRES = true;
if(m_DRTGeometryMetaInfo == NULL) {
// todo
@ -309,23 +247,18 @@ void itkImageProcessor::SetCustom_2Dres(double nX1,double nY1,double nX2,double
ImageType3D::SpacingType Spacing;
Spacing [0] = nX1;
Spacing [1] = nY1;
Spacing [1] = 1.;
Spacing [2] = 1.;
m_DRTGeometryMetaInfo->SetDRT1Spacing(Spacing);
Spacing [0] = nX2;
Spacing [1] = nY2;
Spacing [1] = 1.;
Spacing [2] = 1.;
m_DRTGeometryMetaInfo->SetDRT2Spacing(Spacing);
//TODO UPDATE TO FOLLOW
}
void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2)
{
// image1Size[0] = nX1;
// image1Size[1] = nY1;
// image2Size[0] = nX2;
// image2Size[1] = nY2;
customized_2DSIZE = true;
if(m_DRTGeometryMetaInfo == NULL) {
// todo
@ -334,12 +267,12 @@ void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2)
ImageType3D::SizeType Size;
Size [0] = nX1;
Size [1] = nY1;
Size [1] = 1.;
Size [2] = 1.;
m_DRTGeometryMetaInfo->SetDRT1Size(Size);
Size [0] = nX2;
Size [1] = nY2;
Size [1] = 1.;
Size [2] = 1.;
m_DRTGeometryMetaInfo->SetDRT2Size(Size);
//TODO UPDATE TO FOLLOW
@ -349,8 +282,6 @@ void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2)
int itkImageProcessor::load3DSerie(const char * pcDirName)
{
// std::cout<<"itkImageProcessor::load3DSerie"<<std::endl;
using NamesGeneratorType = itk::GDCMSeriesFileNames;
NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();
@ -412,7 +343,7 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
using ImageIOType = itk::GDCMImageIO;
ImageIOType::Pointer dicomIO = ImageIOType::New();
ImageSeriesReaderType3D::Pointer
imageSeriesReader3D = ImageSeriesReaderType3D::New();
imageSeriesReader3D = ImageSeriesReaderType3D::New();
imageSeriesReader3D->SetImageIO(dicomIO);
imageSeriesReader3D->SetFileNames(fileNames);
@ -440,7 +371,7 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::HFS])){
m_PatOrientation = eImageOrientationType::HFS;
} else if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::FFS])){
m_PatOrientation = eImageOrientationType::FFS;
m_PatOrientation = eImageOrientationType::FFS;
} else {
m_PatOrientation = eImageOrientationType::NotDefined;
}
@ -460,6 +391,8 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
m_VolumeSourceDupli->SetInputImage(caster3D->GetOutput());
m_VolumeSourceDupli->Update();
caster3D = NULL;
}
}
@ -507,10 +440,10 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
/* initialise DRT meta */
m_DRTImage1MetaInfo = DRTImageMetaInformation::New();
m_DRTImage1MetaInfo->SetProjectionAngleLPS(
CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle1IEC())
);
CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle1IEC())
);
m_DRTImage1MetaInfo->SetSCD(
m_DRTGeometryMetaInfo->GetSCD());
@ -519,17 +452,32 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
m_DRTImage2MetaInfo = DRTImageMetaInformation::New();
m_DRTImage2MetaInfo->SetProjectionAngleLPS(
CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle2IEC())
);
CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle2IEC())
);
m_DRTImage2MetaInfo->SetSCD(
m_DRTGeometryMetaInfo->GetSCD());
std::cout<< m_DRTGeometryMetaInfo->GetDRT1Size() <<std::endl;
std::cout<< m_DRTGeometryMetaInfo->GetDRT1Spacing() <<std::endl;
this->UpdateProjectionGeometryMeta();
std::cout<< m_DRTGeometryMetaInfo->GetDRT1Size() <<std::endl;
std::cout<< m_DRTGeometryMetaInfo->GetDRT1Spacing() <<std::endl;
std::cout<< m_DRTImage1MetaInfo->GetSize() <<std::endl;
std::cout<< m_DRTImage1MetaInfo->GetSpacing() <<std::endl;
std::cout<< m_DRTImage1MetaInfo->GetOrigin() <<std::endl;
std::cout<< m_DRTImage2MetaInfo->GetSize() <<std::endl;
std::cout<< m_DRTImage2MetaInfo->GetSpacing() <<std::endl;
std::cout<< m_DRTImage2MetaInfo->GetOrigin() <<std::endl;
// 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
@ -551,38 +499,29 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
m_3DInputChangeInformationToZero->GetOutput();
/* TODO: Here calculate COV
* m_3DProjectionOriginLPSZero
* m_3DProjectionOriginLPS */
return EXIT_SUCCESS;
}
itkImageProcessor::ImageType3D::PointType
itkImageProcessor::CalculateDRTOrigin(
ImageType3D::SizeType m_ImageSize,
ImageType3D::SpacingType m_ImageResolution,
double dSCD
ImageType3D::SizeType m_ImageSize,
ImageType3D::SpacingType m_ImageResolution,
double dSCD
){
ImageType3D::PointType Origin;
double o2Dx, o2Dy;
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 ;
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;
return
Origin;
}
itkImageProcessor::ImageType3D::PointType
@ -614,135 +553,49 @@ itkImageProcessor::CalculateProjectionCenterLPS(
if(bZero)
{
return
m_ProjectionOriginLPS - Origin;
m_ProjectionOriginLPS - Origin;
} else {
return
m_ProjectionOriginLPS;
m_ProjectionOriginLPS;
}
}
//void itkImageProcessor::ApplyVolumeImportTransform(){
// //ImageType3D::PointType ImportOffset;
// if(customized_ImportTransform &&
// customized_RTCouchSetup &&
// customized_RTIsocenter ){
// ImportOffset=
// this->CalcImportVolumeOffset(rtCouchOffset,
// m_CTMetaInfo->GetLPS2IECDirections(),
// rtIsocenterLPS,
// IEC2DCMMapT);
// std::cout<<"ImportOffset> "<<ImportOffset<<std::endl;
// } else {
// ImportOffset.Fill(0.);
// }
// ImageType3D::PointType Origin;
// Origin = m_VolumeSourceDupli->GetOutput()->GetOrigin();
// Origin = Origin - ImportOffset ;
// ChangeInformationFilterType::Pointer
// m_3DInputChangeInformation =
// ChangeInformationFilterType::New();
// m_3DInputChangeInformation->SetInput(
// m_VolumeSourceDupli->GetOutput());
// m_3DInputChangeInformation->SetOutputOrigin(Origin ) ;//NewOrigin);
// m_3DInputChangeInformation->ChangeOriginOn();
// m_3DInputChangeInformation->UpdateOutputInformation();
// m_3DInputChangeInformation->Update();
// /* End of origin change */
// InternalImageType::Pointer m_Image =
// m_3DInputChangeInformation->GetOutput();
// // TODO: move this view origin to new meta variable.
// m_CTMetaInfo->SetOriginLPS(m_Image->GetOrigin());
// /* Calculate center of projection in patient coordinates */
// InternalImageType::DirectionType IECtoLPS_Directions;
// IECtoLPS_Directions =
// m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
// // LPStoIEC_Directions.GetTranspose();
// if(customized_ProjectionCenter)
// {
// m_3DProjectionOriginLPS =
// IECtoLPS_Directions * m_DRTGeometryMetaInfo->GetProjectionCenter();
// /* in longitudinal direction (Sup-Inf), we put it in the
// * middle of the volume */
// m_3DProjectionOriginLPS[2] = m_CTMetaInfo->GetCOV()[2];
// } else {
// m_3DProjectionOriginLPS = m_CTMetaInfo->GetCOV();
// }
// /* calculate also the same center with zero origin
// * This is required to bring back the DRT into LPS */
// m_3DProjectionOriginLPSZero =
// m_3DProjectionOriginLPS - m_CTMetaInfo->GetOriginLPS() ;
// /* We should know everything we need to calculate the
// * origin and direction of projection images */
// //std::cout<< "m_3DProjectionOriginLPS "<<m_3DProjectionOriginLPS <<std::endl;
// //std::cout<< "m_3DProjectionOriginLPSZero "<<m_3DProjectionOriginLPSZero <<std::endl;
// if(true){
// std::cout<<" -------- 3D IMAGE --------"<<std::endl;
//// std::cout<<"size3D " <<size3D <<std::endl;
//// std::cout<<"resolution3D " <<resolution3D <<std::endl;
//// std::cout<<"imagDirection " <<imagDirection <<std::endl;
//// std::cout<<"First Voxel position (origin) " <<m_3DOriginLPS <<std::endl;
//// std::cout<<"Last Voxel position: "<<LastVoxelPosition<<std::endl;
// std::cout<<"Image 3D COV: "<< m_CTMetaInfo->GetCOV() <<std::endl;
// std::cout<<"Data collection center LPS: "<<m_3DProjectionOriginLPS <<std::endl;
// std::cout<<"Data collection center LPS (Zero Origin): "<<m_3DProjectionOriginLPSZero<<std::endl;
// std::cout<<" -------- -------- --------"<<std::endl;
// /*std::cout<<"META Proj Origin LPS"<<
// m_DRTImage1MetaInfo->GetProjectionOriginLPS()<<std::endl;
// std::cout<<"META Proj Origin LPS"<<
// m_DRTImage1MetaInfo->GetProjectionOriginLPSZero()<<std::endl;
// std::cout<<"META Proj Origin LPS"<<
// m_DRTImage2MetaInfo->GetProjectionOriginLPS()<<std::endl;
// std::cout<<"META Proj Origin LPS"<<
// m_DRTImage2MetaInfo->GetProjectionOriginLPSZero()<<std::endl;*/
// }
//}
const double* itkImageProcessor::GetRTImportOffset()
const std::vector <double> itkImageProcessor::GetRTImportOffset()
{
std::vector <double> vOffset;
vOffset.clear();
vOffset.push_back(
m_CTMetaInfo->GetImportOffset()[0]
);
vOffset.push_back(
m_CTMetaInfo->GetImportOffset()[1]
);
vOffset.push_back(
m_CTMetaInfo->GetImportOffset()[2]
);
return
ImportOffset.GetDataPointer();
vOffset;
}
const std::vector <double> itkImageProcessor::GetLPStoProjectionGeoLPSOffset()
{
std::vector <double> vOffset;
vOffset.clear();
vOffset.push_back(
m_DRTImage1MetaInfo->GetLPStoProjectionGeoLPSOffset()[0]
);
);
vOffset.push_back(
m_DRTImage1MetaInfo->GetLPStoProjectionGeoLPSOffset()[1]
);
);
vOffset.push_back(
m_DRTImage1MetaInfo->GetLPStoProjectionGeoLPSOffset()[2]
);
);
return
vOffset;
}
double
@ -1052,57 +905,57 @@ double itkImageProcessor::GetLocalizerDisplayWindowLevel(int iImg)
switch (iImg) {
case 0:
case 0:
m_TImageMeta = m_TImage1MetaInfo;
break;
break;
case 1:
case 1:
m_TImageMeta = m_TImage2MetaInfo;
break;
break;
default:
return 0;
break;
default:
return 0;
break;
}
}
if(m_TImageMeta == NULL){
return 0;
} else {
return
m_TImageMeta->GetWLLevel();
}
if(m_TImageMeta == NULL){
return 0;
} else {
return
m_TImageMeta->GetWLLevel();
}
}
double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg)
double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg)
{
TopogramImageMetaInformation::Pointer m_TImageMeta;
switch (iImg) {
switch (iImg) {
case 0:
m_TImageMeta = m_TImage1MetaInfo;
break;
case 0:
m_TImageMeta = m_TImage1MetaInfo;
break;
case 1:
m_TImageMeta = m_TImage2MetaInfo;
case 1:
m_TImageMeta = m_TImage2MetaInfo;
break;
break;
default:
return 0;
break;
default:
return 0;
break;
}
}
if(m_TImageMeta == NULL){
return 0;
} else {
return
m_TImageMeta->GetWLWindow();
}
if(m_TImageMeta == NULL){
return 0;
} else {
return
m_TImageMeta->GetWLWindow();
}
}
@ -1122,13 +975,13 @@ void itkImageProcessor::InitializeProjector()
const double dtr = (atan(1.0) * 4.0) / 180.0;
transform->SetRotation(
dtr * RZero[0],
dtr * RZero[1],
dtr * RZero[2]);
dtr * RZero[1],
dtr * RZero[2]);
if(m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() !=
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ) {
// TODO
}
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ) {
// TODO
}
transform->SetCenter(
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() );
@ -1142,9 +995,7 @@ void itkImageProcessor::InitializeProjector()
interpolator1->SetThreshold(
m_DRTGeometryMetaInfo->GetIntensityThreshold()
);
interpolator1->SetTransform(transform);
interpolator1->Initialize();
// 2D Image 2
@ -1157,7 +1008,6 @@ void itkImageProcessor::InitializeProjector()
m_DRTGeometryMetaInfo->GetIntensityThreshold()
);
interpolator2->SetTransform(transform);
interpolator2->Initialize();
}
@ -1229,9 +1079,9 @@ void itkImageProcessor::loadRTPlanInfo(
void itkImageProcessor::UpdateProjectionGeometryMeta(){
if(m_CTMetaInfo == NULL){
//TODO.
}
if(m_CTMetaInfo == NULL){
//TODO.
}
m_DRTImage1MetaInfo->SetProjectionOriginLPS(
CalculateProjectionCenterLPS(
@ -1339,59 +1189,15 @@ void itkImageProcessor::UpdateProjectionGeometryMeta(){
To be very careful editing ...
*/
itkImageProcessor::ImageType3D::PointType
itkImageProcessor::CalcDRTImageOrigin(itkImageProcessor::InternalImageType::Pointer m_Image,
itkImageProcessor::ImageType3D::PointType m_VolCOOV,
double dAngle,
InternalImageType::DirectionType stdDRT2LPS
){
using ImageRegionType3D = itkImageProcessor::ImageType3D::RegionType;
using SizeType3D = ImageRegionType3D::SizeType;
using ImageDirectionType3D = itkImageProcessor::ImageType3D::DirectionType;
using VectorType = itk::Vector<double, 3>;
ImageRegionType3D region3D;
SizeType3D size3D;
itk::Vector<double, 3> 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, stdDRT2LPS);
ImageType3D::PointType NewOrigin =
m_VolCOOV + DRT2LPS * (m_Image ->GetOrigin() - DRTCOV);
return
NewOrigin;
}
itkImageProcessor::ImageType3D::PointType
itkImageProcessor::CalcDRTImageOrigin(
ImageType3D::PointType m_DRTOrigin,
ImageType3D::PointType m_DRTCOV,
ImageType3D::PointType m_ProjectionCenterLPS,
double dAngle,
InternalImageType::DirectionType stdDRT2LPS
){
ImageType3D::PointType m_DRTOrigin,
ImageType3D::PointType m_DRTCOV,
ImageType3D::PointType m_ProjectionCenterLPS,
double dAngle,
InternalImageType::DirectionType stdDRT2LPS
){
itkImageProcessor::InternalImageType::DirectionType DRT2LPS
= this->CalcDRTImageDirections(dAngle, stdDRT2LPS);
@ -1423,8 +1229,8 @@ void itkImageProcessor::GetProjectionImages(){
// The transform is determined by the parameters and the rotation center.
if(m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() !=
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ) {
// TODO
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() ) {
// TODO
}
finalTransform->SetCenter(
@ -1448,6 +1254,7 @@ void itkImageProcessor::GetProjectionImages(){
interpolator1->Initialize();
resampleFilter1->SetInterpolator(interpolator1);
resampleFilter1->SetSize(m_DRTImage1MetaInfo->GetSize());
resampleFilter1->SetOutputSpacing(m_DRTImage1MetaInfo->GetSpacing());
resampleFilter1->SetOutputOrigin(m_DRTImage1MetaInfo->GetOrigin());
@ -1471,11 +1278,9 @@ void itkImageProcessor::GetProjectionImages(){
resampleFilter2->SetOutputSpacing(m_DRTImage2MetaInfo->GetSpacing());
resampleFilter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOrigin());
resampleFilter1->Update();
resampleFilter2->Update();
/* here probably we should move everything back to dicom */
/* Projection images are in a fancy reference system...
@ -1511,6 +1316,7 @@ void itkImageProcessor::GetProjectionImages(){
filter2->ChangeOriginOn();
filter2->UpdateOutputInformation();
filter1->Update();
filter2->Update();
@ -1552,60 +1358,60 @@ itkImageProcessor::CalcImportVolumeOffset(
void itkImageProcessor::Write2DImages(){
using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
using WriterType = itk::ImageFileWriter<OutputImageType>;
// using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
// using WriterType = itk::ImageFileWriter<OutputImageType>;
if(image2D1Loaded)
{
// if(image2D1Loaded)
// {
// Rescale the intensity of the projection images to 0-255 for output.
RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
rescaler1->SetOutputMinimum(0);
rescaler1->SetOutputMaximum(255);
rescaler1->SetInput(m_PASourceDupli->GetOutput());
// // Rescale the intensity of the projection images to 0-255 for output.
// RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
// rescaler1->SetOutputMinimum(0);
// rescaler1->SetOutputMaximum(255);
// rescaler1->SetInput(m_PASourceDupli->GetOutput());
WriterType::Pointer writer1 = WriterType::New();
// WriterType::Pointer writer1 = WriterType::New();
writer1->SetFileName("/mnt/gdrive/Scratch/gfattori/2D1.tif");
writer1->SetInput(rescaler1->GetOutput());
// writer1->SetFileName("/mnt/gdrive/Scratch/gfattori/2D1.tif");
// writer1->SetInput(rescaler1->GetOutput());
try
{
std::cout << "Writing image 2D1" << std::endl;
writer1->Update();
}
catch (itk::ExceptionObject & err)
{
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
// try
// {
// std::cout << "Writing image 2D1" << std::endl;
// writer1->Update();
// }
// catch (itk::ExceptionObject & err)
// {
// std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
// std::cerr << err << std::endl;
// }
}
// }
if(image2D2Loaded){
// Rescale the intensity of the projection images to 0-255 for output.
RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
rescaler2->SetOutputMinimum(0);
rescaler2->SetOutputMaximum(255);
rescaler2->SetInput(m_LATSourceDupli->GetOutput());
// if(image2D2Loaded){
// // Rescale the intensity of the projection images to 0-255 for output.
// RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
// rescaler2->SetOutputMinimum(0);
// rescaler2->SetOutputMaximum(255);
// rescaler2->SetInput(m_LATSourceDupli->GetOutput());
WriterType::Pointer writer2 = WriterType::New();
// WriterType::Pointer writer2 = WriterType::New();
writer2->SetFileName("/mnt/gdrive/Scratch/gfattori/2D2.tif");
writer2->SetInput(rescaler2->GetOutput());
// writer2->SetFileName("/mnt/gdrive/Scratch/gfattori/2D2.tif");
// writer2->SetInput(rescaler2->GetOutput());
try
{
std::cout << "Writing image 2D2" << std::endl;
writer2->Update();
}
catch (itk::ExceptionObject & err)
{
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
}
// try
// {
// std::cout << "Writing image 2D2" << std::endl;
// writer2->Update();
// }
// catch (itk::ExceptionObject & err)
// {
// std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
// std::cerr << err << std::endl;
// }
// }
}
@ -1788,36 +1594,31 @@ vtkImageData* itkImageProcessor::GetProjection1VTK()
vtkMatrix4x4 * itkImageProcessor::GetProjection1VTKTransform()
{
// std::cout<< "Composed " <<std::endl;
// interpolator1->GetComposedTransform()->Print(std::cout);
// std::cout<< "Inverse " <<std::endl;
// interpolator1->GetInverseTransform()->Print(std::cout);
// std::cout<< "Composed " <<std::endl;
// interpolator1->GetComposedTransform()->Print(std::cout);
// std::cout<< "Inverse " <<std::endl;
// interpolator1->GetInverseTransform()->Print(std::cout);
m_Projection1VTKTransform->Identity();
for(int iIdx = 0; iIdx<3 ; iIdx++){
for(int jIdx = 0; jIdx<3 ; jIdx++){
m_Projection1VTKTransform->SetElement(iIdx,jIdx,
interpolator1->GetInverseTransform()->GetMatrix().GetVnlMatrix()[iIdx][jIdx]);
interpolator1->GetInverseTransform()->GetMatrix().GetVnlMatrix()[iIdx][jIdx]);
}
}
m_Projection1VTKTransform->SetElement(0,3,
interpolator1->GetInverseTransform()->GetTranslation()[0]);
interpolator1->GetInverseTransform()->GetTranslation()[0]);
m_Projection1VTKTransform->SetElement(1,3,
interpolator1->GetInverseTransform()->GetTranslation()[1]);
interpolator1->GetInverseTransform()->GetTranslation()[1]);
m_Projection1VTKTransform->SetElement(2,3,
interpolator1->GetInverseTransform()->GetTranslation()[2]);
//m_Projection1VTKTransform->Print(std::cout);
interpolator1->GetInverseTransform()->GetTranslation()[2]);
return
m_Projection1VTKTransform;
}
vtkMatrix4x4 * itkImageProcessor::GetProjection2VTKTransform()
@ -1828,18 +1629,18 @@ vtkMatrix4x4 * itkImageProcessor::GetProjection2VTKTransform()
for(int iIdx = 0; iIdx<3 ; iIdx++){
for(int jIdx = 0; jIdx<3 ; jIdx++){
m_Projection2VTKTransform->SetElement(iIdx,jIdx,
interpolator2->GetInverseTransform()->GetMatrix().GetVnlMatrix()[iIdx][jIdx]);
interpolator2->GetInverseTransform()->GetMatrix().GetVnlMatrix()[iIdx][jIdx]);
}
}
m_Projection2VTKTransform->SetElement(0,3,
interpolator2->GetInverseTransform()->GetTranslation()[0]);
interpolator2->GetInverseTransform()->GetTranslation()[0]);
m_Projection2VTKTransform->SetElement(1,3,
interpolator2->GetInverseTransform()->GetTranslation()[1]);
interpolator2->GetInverseTransform()->GetTranslation()[1]);
m_Projection2VTKTransform->SetElement(2,3,
interpolator2->GetInverseTransform()->GetTranslation()[2]);
interpolator2->GetInverseTransform()->GetTranslation()[2]);
@ -1977,7 +1778,6 @@ void itkImageProcessor::SetInitialTranslations(double dX, double dY, double dZ)
IECtoLPS_Directions =
m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
// LPStoIEC_Directions.GetTranspose();
ImageType3D::PointType LPSTranslations =
IECtoLPS_Directions * IECTranslations;
@ -1986,9 +1786,6 @@ void itkImageProcessor::SetInitialTranslations(double dX, double dY, double dZ)
TZero[1]= LPSTranslations[1];
TZero[2]= LPSTranslations[2];
// std::cout<< "itkImageProcessor::SetInitialTranslations IEC : "<<IECTranslations <<std::endl;
// std::cout<< "itkImageProcessor::SetInitialTranslations LPS : "<<TZero[0]<<" "<<TZero[1]<<" "<<TZero[2] <<std::endl;
}
void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ)

63
reg23Topograms/itkDTRrecon/itkImageProcessor.h
View File

@ -97,8 +97,6 @@ public:
/** Custom settings of the projection images */
void SetCustom_2Dres(double,double,double,double);
void SetCustom_2Dsize(int,int,int,int);
//void SetCustom_RTIsocenter_LPS(double, double, double);
//void SetCustom_RTCouchSetup_IEC(double, double, double);
void SetCustom_ImportTransform(double, double, double,
double, double, double);
void SetCustom_ProjectionCenterFixedAxes_IEC(double,double,double);
@ -114,7 +112,7 @@ public:
const std::vector <double> GetLPStoProjectionGeoLPSOffset();
/** Get Projection origin in LPS coordinates*/
const double* GetRTImportOffset();
const std::vector <double> GetRTImportOffset();
/** Apply transform to CT volume.
@ -211,11 +209,6 @@ private:
ImageType3D::PointType m_ProjectionCenter,
bool bZero);
ImageType3D::PointType ImportOffset;
ImageType3D :: PointType
ConversionOffset;
TransformType::Pointer
transform;
@ -223,7 +216,6 @@ private:
interpolator1,
interpolator2;
ITKtoVTKFilterType::Pointer
toVTK2D1,
toVTK2D2,
@ -243,10 +235,6 @@ private:
ChangeInformationFilterType::Pointer
m_3DInputChangeInformationToZero;
// ImageType3D :: PointType
// m_3DProjectionOriginLPS,
// m_3DProjectionOriginLPSZero;
void UpdateProjectionGeometryMeta();
@ -257,13 +245,6 @@ private:
CalcDRTImageDirections(double angle,
InternalImageType::DirectionType DRT2LPS);
ImageType3D::PointType
CalcDRTImageOrigin(InternalImageType::Pointer m_Image,
ImageType3D::PointType m_VolCOOV,
double dAngle,
InternalImageType::DirectionType stdDRT2LPS
);
ImageType3D::PointType
CalcDRTImageOrigin(
ImageType3D::PointType m_DRTOrigin,
@ -271,7 +252,7 @@ private:
ImageType3D::PointType m_ProjectionCenterLPS,
double dAngle,
InternalImageType::DirectionType stdDRT2LPS
);
);
double
@ -286,49 +267,12 @@ private:
ImageType3D::PointType rtIsocenterLPS,
ImageType3D::PointType IEC2DCMMapT);
// double image1res[2], image2res[2];
// int image1Size[2], image2Size[2];
double TZero[3], RZero[3];
bool
customized_2DRES,
customized_2DSIZE,
customized_RTIsocenter,
customized_RTCouchSetup,
customized_ImportTransform,
customized_ProjectionCenter;
bool
image2D1Loaded,
image2D2Loaded,
image3DLoaded;
/* Transform between IEC Support and
* IEC Scanner frame of reference */
// ImageType3D::PointType
// IEC2DCMMapT,
// IEC2DCMMapR;
/* RT Plan isocenter in LPS coordinates */
// ImageType3D::PointType
// rtIsocenterLPS;
/* RT Plan couch SETUP offset corresponding to
* the RT Plan isocenter */
// ImageType3D::PointType
// rtCouchOffset;
InternalImageType::DirectionType
Std_DRT2LPS;
vtkMatrix4x4
* m_Projection1VTKTransform,
* m_Projection2VTKTransform;
@ -355,6 +299,9 @@ private:
RTGeometryMetaInformation::Pointer
m_RTMetaInfo;
TransformMetaInformation::Pointer
m_TransformMetaInfo;
};