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

Initial code cleanup. Removing AutoReg functions and better split of itk pipeline for calculation of LPS coordinates and actual projection

Browse Source
This commit is contained in:
fattori_g
2022-02-11 23:54:42 +01:00
parent 70498e03b6
commit a65834c60e
2 changed files with 287 additions and 314 deletions
Download Patch File Download Diff File Expand all files Collapse all files

427
reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
View File

@ -29,13 +29,6 @@ gfattori 08.11.2021
/* Settings */
#define CT_IMPORTOFFSET_TXIEC 0
#define CT_IMPORTOFFSET_TYIEC 0
#define CT_IMPORTOFFSET_TZIEC 0
#define CT_IMPORTOFFSET_RXIEC 0
#define CT_IMPORTOFFSET_RYIEC 0
#define CT_IMPORTOFFSET_RZIEC 0
#define CT_DEFAULT_TABLEHEIGHT 175
/* From Files */
@ -46,9 +39,6 @@ gfattori 08.11.2021
#define RT_LNG 0
#define RT_VRT 0
#define TMP_CT3X 4.1305
#define TMP_CT3Y -0.9476
#define TMP_CT3Z 252.7127
namespace itk
{
@ -114,29 +104,27 @@ itkImageProcessor::itkImageProcessor()
interpolator1 = InterpolatorType::New();
interpolator2 = InterpolatorType::New();
imageReader3D = ImageReaderType3D::New();
imageSeriesReader3D = ImageSeriesReaderType3D::New();
//imageSeriesReader3D = ImageSeriesReaderType3D::New();
scd = 1000.; // Source to isocenter distance
threshold = 0.;
d_scd = 1000.; // Source to isocenter distance
d_intensityThreshold = 0.;
projAngle1 = -999.;
projAngle2 = -999.;
isocenter[0]=isocenter[1]=isocenter[2] = 0.;
image1res[0]=image2res[1]= 2.;
image1center[0]=image1center[1]= 0.;
image2center[0]=image2center[1]= 0.;
//image1center[0]=image1center[1]= 0.;
//image2center[0]=image2center[1]= 0.;
image1Size[0] = image1Size[1] = 512;
image2Size[0] = image2Size[1] = 512;
// image3DCOV[0] = image3DCOV[1] = image3DCOV[2] = 0.;
image1IntensityWindow[0] = image1IntensityWindow[1] = 0.;
image2IntensityWindow[0] = image2IntensityWindow[1] = 0.;
d_image1IntensityWindow[0] = d_image1IntensityWindow[1] = 0.;
d_image2IntensityWindow[0] = d_image2IntensityWindow[1] = 0.;
TZero[0]=TZero[1]=TZero[2]=0.;
RZero[0]=RZero[1]=RZero[2]=0.;
customized_iso = false;
customized_2DCX = false;
// customized_2DCX = false;
customized_2DRES = false;
customized_2DSIZE = false;
@ -144,9 +132,8 @@ itkImageProcessor::itkImageProcessor()
image2D2Loaded= false;
image3DLoaded= false;
caster3D = CastFilterType3D::New();
caster2D1 = CastFilterType3D::New();
caster2D2 = CastFilterType3D::New();
// caster3D = CastFilterType3D::New();
toVTK2D1 = ITKtoVTKFilterType::New();
toVTK2D2 = ITKtoVTKFilterType::New();
@ -157,7 +144,7 @@ itkImageProcessor::itkImageProcessor()
m_LATSourceDupli = DuplicatorType::New();
m_PASourceDupli = DuplicatorType::New();
m_VolumeSourceDupli = DuplicatorType::New();
for(int iIdx = 0 ; iIdx < 3; iIdx++){
for(int jIdx = 0 ; jIdx < 3; jIdx++){
@ -167,14 +154,8 @@ itkImageProcessor::itkImageProcessor()
m_3DInputChangeInformation = ChangeInformationFilterInputType::New();
d_CTImportOffset_txiec = CT_IMPORTOFFSET_TXIEC;
d_CTImportOffset_tyiec = CT_IMPORTOFFSET_TYIEC;
d_CTImportOffset_tziec = CT_IMPORTOFFSET_TZIEC;
d_CTImportOffset_rxiec = CT_IMPORTOFFSET_RXIEC;
d_CTImportOffset_ryiec = CT_IMPORTOFFSET_RYIEC;
d_CTImportOffset_rziec = CT_IMPORTOFFSET_RZIEC;
m_3DInputChangeInformationToZero =
ChangeInformationFilterType::New();
}
@ -191,12 +172,12 @@ void itkImageProcessor::PrintSelf(std::ostream& os, Indent indent) const
void itkImageProcessor::SetIntensityThreshold(double dT)
{
threshold = dT;
d_intensityThreshold = dT;
}
void itkImageProcessor::SetSCD(double dDist)
{
scd = dDist;
d_scd = dDist;
}
void itkImageProcessor::SetCustom_Isocenter(double dX,double dY,double dZ)
@ -225,14 +206,14 @@ void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2)
customized_2DSIZE = true;
}
void itkImageProcessor::SetCustom_2Dcenter(double dX1,double dY1, double dX2,double dY2)
{
image1center[0] = dX1;
image1center[1] = dY1;
image2center[0] = dX2;
image2center[1] = dY2;
customized_2DCX = true;
}
//void itkImageProcessor::SetCustom_2Dcenter(double dX1,double dY1, double dX2,double dY2)
//{
// image1center[0] = dX1;
// image1center[1] = dY1;
// image2center[0] = dX2;
// image2center[1] = dY2;
// customized_2DCX = true;
//}
int itkImageProcessor::load3DSerie(const char * pcDirName)
{
@ -295,6 +276,9 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
using ImageIOType = itk::GDCMImageIO;
ImageIOType::Pointer dicomIO = ImageIOType::New();
ImageSeriesReaderType3D::Pointer
imageSeriesReader3D = ImageSeriesReaderType3D::New();
imageSeriesReader3D->SetImageIO(dicomIO);
imageSeriesReader3D->SetFileNames(fileNames);
imageSeriesReader3D->SetNumberOfWorkUnits(8);
@ -302,9 +286,6 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
imageSeriesReader3D->Update();
// eImageOrientationType
// currImgOrient;
/* check patient orientation */
if(fileNames.size()>0){
gdcm::Reader R;
@ -322,19 +303,20 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
*std::remove(sTmpString, sTmpString + strlen(sTmpString), ' ') = 0;
if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::HFS])){
std::cout<<"Image Orientation: HFS"<<std::endl;
// currImgOrient = eImageOrientationType::HFS;
//std::cout<<"Image Orientation: HFS"<<std::endl;
m_3DPatOrientation = eImageOrientationType::HFS;
} else if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::FFS])){
std::cout<<"Image Orientation: FFS"<<std::endl;
// currImgOrient = eImageOrientationType::FFS;
//std::cout<<"Image Orientation: FFS"<<std::endl;
m_3DPatOrientation = eImageOrientationType::FFS;
} else {
m_3DPatOrientation = eImageOrientationType::NotDefined;
std::cerr<< "Image Orientation: Unrecognised"<< sTmpString <<std::endl;
}
free (sTmpString);
} else {
return -1;
}
/* save the current LPS-to-IEC */
@ -362,156 +344,41 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
projAngle2 = CalcProjectionAngleLPS(m_3DPatOrientation, projAngle2_IEC);
CastFilterType3D::Pointer caster3D =
CastFilterType3D::New();
caster3D->SetInput(imageSeriesReader3D->GetOutput());
caster3D->Update();
/* TODO: here is the right place to change origin when required
* after this point, reference info are saved */
ImageType3D::PointType rtCouchOffset;
rtCouchOffset[0]=9.60802;
rtCouchOffset[1]=1333.925;
rtCouchOffset[2]=12.2;
ImageType3D::PointType rtIsocenterLPS;
rtIsocenterLPS[0]= 9.01502;
rtIsocenterLPS[1]= -157.708;
rtIsocenterLPS[2]= -1186.84;
ImageType3D::PointType IEC2DCMMapT;
IEC2DCMMapT[0]= 4.10766744;
IEC2DCMMapT[1]= -2316.44798289;
IEC2DCMMapT[2]= 172.82908123;
ImageType3D::PointType ImportOffset=
this->CalcImportVolumeOffset(rtCouchOffset,
LPStoIEC_Directions,
rtIsocenterLPS,
IEC2DCMMapT);
ImageType3D::PointType Origin;
Origin = imageSeriesReader3D->GetOutput()->GetOrigin();
std::cout<< "OldOrigin: " <<Origin<<std::endl;
Origin = Origin - ImportOffset ;
std::cout<< "NewOrigin: " <<Origin<<std::endl;
// exit(1);
// ImageType3D::PointType NewOrigin;
// ImageType3D::PointType Offset;
// Offset[0] = -TMP_CT3X;
// Offset[1] = -TMP_CT3Y;
// Offset[2] = -TMP_CT3Z;
// NewOrigin = imageSeriesReader3D->GetOutput()->GetOrigin();
// std::cout<< "OldOrigin: " <<NewOrigin<<std::endl;
// NewOrigin = NewOrigin - Offset;
m_3DInputChangeInformation->SetInput(imageSeriesReader3D->GetOutput());
m_3DInputChangeInformation->SetOutputOrigin(Origin ) ;//NewOrigin);
m_3DInputChangeInformation->ChangeOriginOn();
m_3DInputChangeInformation->UpdateOutputInformation();
m_3DInputChangeInformation->Update();
std::cout<< "NewOrigin: " <<m_3DInputChangeInformation->GetOutput()->GetOrigin() <<std::endl;
/* End of origin change */
image3DIn = m_3DInputChangeInformation->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<double, 3> resolution3D = image3DIn->GetSpacing();
ImageDirectionType3D imagDirection = image3DIn->GetDirection();
/* 3D volume origin in patient coordinates */
m_3DOriginLPS = image3DIn->GetOrigin();
/* calculate image size in 3D Space by finding the last voxel position */
using VectorType = itk::Vector<double, 3>;
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 =
m_3DOriginLPS + (imagDirection * Dest);
ImageType3D::PointType image3DCOVLPS;
image3DCOVLPS[0] = (m_3DOriginLPS[0]+LastVoxelPosition[0])/2;
image3DCOVLPS[1] = (m_3DOriginLPS[1]+LastVoxelPosition[1])/2;
image3DCOVLPS[2] = (m_3DOriginLPS[2]+LastVoxelPosition[2])/2;
/* Calculate center of projection in patient coordinates */
/* center of projection in an IEC reference at Patient Origin */
ImageType3D::PointType ProjectionCenterFixedAxes;
ProjectionCenterFixedAxes[0] = 0.;
ProjectionCenterFixedAxes[1] = 0.;
ProjectionCenterFixedAxes[2] = CT_DEFAULT_TABLEHEIGHT;
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions = LPStoIEC_Directions.GetTranspose();
m_3DProjectionOriginLPS =
IECtoLPS_Directions * ProjectionCenterFixedAxes ;
/* in longitudinal direction (Sup-Inf), we put it in the
* middle of the volume */
m_3DProjectionOriginLPS[2] = image3DCOVLPS[2];
/* calculate also the same center with zero origin
* This is required to bring back the DRT into LPS */
m_3DProjectionOriginLPSZero =
m_3DProjectionOriginLPS - m_3DOriginLPS;
/* We should know everything we need to calculate the
* origin and direction of projection images */
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: "<<image3DCOVLPS <<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;
}
m_VolumeSourceDupli->SetInputImage(caster3D->GetOutput());
m_VolumeSourceDupli->Update();
// To simply Siddon-Jacob's fast ray-tracing algorithm, we force the origin of the CT image
// to be (0,0,0). Because we align the CT isocenter with the central axis, the projection
// geometry is fully defined. The origin of the CT image becomes irrelavent.
ImageType3D::PointType image3DOrigin;
image3DOrigin[0] = 0.0;
image3DOrigin[1] = 0.0;
image3DOrigin[2] = 0.0;
image3DIn->SetOrigin(image3DOrigin);
ImageType3D::PointType pZeroOrigin;
pZeroOrigin[0] = 0.;
pZeroOrigin[1] = 0.;
pZeroOrigin[2] = 0.;
m_3DInputChangeInformationToZero->SetInput(
m_VolumeSourceDupli->GetOutput());
m_3DInputChangeInformationToZero->SetOutputOrigin(
pZeroOrigin);
m_3DInputChangeInformationToZero->UpdateOutputInformation();
m_3DInputChangeInformationToZero->Update();
image3DIn =
m_3DInputChangeInformationToZero->GetOutput();
caster3D->SetInput(image3DIn);
caster3D->Update();
/* TODO: Here calculate COV
* m_3DProjectionOriginLPSZero
* m_3DProjectionOriginLPS */
}
}
catch (itk::ExceptionObject & ex)
{
@ -527,6 +394,130 @@ int itkImageProcessor::load3DSerie(const char * pcDirName)
}
void itkImageProcessor::ApplyVolumeImportTransform(){
/* TODO: here is the right place to change origin when required
* after this point, reference info are saved */
ImageType3D::PointType rtCouchOffset;
rtCouchOffset[0]=9.60802;
rtCouchOffset[1]=1333.925;
rtCouchOffset[2]=12.2;
ImageType3D::PointType rtIsocenterLPS;
rtIsocenterLPS[0]= 9.01502;
rtIsocenterLPS[1]= -157.708;
rtIsocenterLPS[2]= -1186.84;
ImageType3D::PointType IEC2DCMMapT;
IEC2DCMMapT[0]= 4.10766744;
IEC2DCMMapT[1]= -2316.44798289;
IEC2DCMMapT[2]= 172.82908123;
ImageType3D::PointType ImportOffset=
this->CalcImportVolumeOffset(rtCouchOffset,
LPStoIEC_Directions,
rtIsocenterLPS,
IEC2DCMMapT);
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();
/* 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 = m_Image->GetBufferedRegion();
SizeType3D size3D = region3D.GetSize();
// ImageType3D::PointType origin3D = image3DIn->GetOrigin();
const itk::Vector<double, 3> resolution3D = m_Image ->GetSpacing();
ImageDirectionType3D imagDirection = m_Image ->GetDirection();
/* 3D volume origin in patient coordinates */
m_3DOriginLPS = m_Image ->GetOrigin();
/* calculate image size in 3D Space by finding the last voxel position */
using VectorType = itk::Vector<double, 3>;
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 =
m_3DOriginLPS + (imagDirection * Dest);
ImageType3D::PointType image3DCOVLPS;
image3DCOVLPS[0] = (m_3DOriginLPS[0]+LastVoxelPosition[0])/2;
image3DCOVLPS[1] = (m_3DOriginLPS[1]+LastVoxelPosition[1])/2;
image3DCOVLPS[2] = (m_3DOriginLPS[2]+LastVoxelPosition[2])/2;
/* Calculate center of projection in patient coordinates */
/* center of projection in an IEC reference at Patient Origin */
ImageType3D::PointType ProjectionCenterFixedAxes;
ProjectionCenterFixedAxes[0] = 0.;
ProjectionCenterFixedAxes[1] = 0.;
ProjectionCenterFixedAxes[2] = CT_DEFAULT_TABLEHEIGHT;
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions = LPStoIEC_Directions.GetTranspose();
m_3DProjectionOriginLPS =
IECtoLPS_Directions * ProjectionCenterFixedAxes ;
/* in longitudinal direction (Sup-Inf), we put it in the
* middle of the volume */
m_3DProjectionOriginLPS[2] = image3DCOVLPS[2];
/* calculate also the same center with zero origin
* This is required to bring back the DRT into LPS */
m_3DProjectionOriginLPSZero =
m_3DProjectionOriginLPS - m_3DOriginLPS;
/* We should know everything we need to calculate the
* origin and direction of projection images */
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: "<<image3DCOVLPS <<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;
}
}
double
itkImageProcessor::CalcProjectionAngleLPS(
tPatOrientation pOrient,
@ -771,7 +762,7 @@ int itkImageProcessor::load2D(const char * pcFName){
switch (currProjOrientation) {
case eProjectionOrientationType::PA:
m_ImageIntensity = image1IntensityWindow;
m_ImageIntensity = d_image1IntensityWindow;
m_Duplicator = m_PASourceDupli;
// m_ImLoaded = &image2D1Loaded;
@ -779,7 +770,7 @@ int itkImageProcessor::load2D(const char * pcFName){
case eProjectionOrientationType::LAT:
m_ImageIntensity = image2IntensityWindow;
m_ImageIntensity = d_image2IntensityWindow;
m_Duplicator = m_LATSourceDupli;
// m_ImLoaded = &image2D2Loaded;
@ -803,17 +794,17 @@ int itkImageProcessor::load2D(const char * pcFName){
}
double* itkImageProcessor::GetImageIntensityWindow(int iImg)
const double* itkImageProcessor::GetImageIntensityWindow(int iImg)
{
switch (iImg) {
case 0:
return image1IntensityWindow;
return d_image1IntensityWindow;
break;
case 1:
return image2IntensityWindow;
return d_image2IntensityWindow;
break;
default:
@ -823,12 +814,11 @@ double* itkImageProcessor::GetImageIntensityWindow(int iImg)
}
}
#include "itkLinearInterpolateImageFunction.h"
void itkImageProcessor::InitializeProjector()
{
caster3D->Update();
transform->SetComputeZYX(true);
TransformType::OutputVectorType translation;
@ -838,31 +828,25 @@ void itkImageProcessor::InitializeProjector()
transform->SetTranslation(translation);
const double dtr = (atan(1.0) * 4.0) / 180.0;
transform->SetRotation(dtr * RZero[0], dtr * RZero[1], dtr * RZero[2]);
transform->SetRotation(
dtr * RZero[0],
dtr * RZero[1],
dtr * RZero[2]);
isocenter3D = m_3DProjectionOriginLPSZero;
transform->SetCenter(m_3DProjectionOriginLPSZero);
// std::cout<<" -------- PROJECTOR GEOMETRY --------"<<std::endl;
std::cout<<"Projector Isocenter "<<isocenter3D[0]<<" "<<isocenter3D[1]<<" "<<isocenter3D[2] <<std::endl;
// std::cout<<origin3D[0]<<" "<<resolution3D[0]<<" "<<static_cast<double>(size3D[0])<<" "<<static_cast<double>(size3D[0]) / 2.0<<std::endl;
// std::cout<<"Projector scd "<<scd <<std::endl;
// std::cout<<" -------- -------- --------"<<std::endl;
transform->SetCenter(isocenter3D);
// 2D Image 1
interpolator1->SetProjectionAngle(dtr * projAngle1);
interpolator1->SetFocalPointToIsocenterDistance(scd);
interpolator1->SetThreshold(threshold);
interpolator1->SetFocalPointToIsocenterDistance(d_scd);
interpolator1->SetThreshold(d_intensityThreshold);
interpolator1->SetTransform(transform);
interpolator1->Initialize();
// 2D Image 2
interpolator2->SetProjectionAngle(dtr * projAngle2);
interpolator2->SetFocalPointToIsocenterDistance(scd);
interpolator2->SetThreshold(threshold);
interpolator2->SetFocalPointToIsocenterDistance(d_scd);
interpolator2->SetThreshold(d_intensityThreshold);
interpolator2->SetTransform(transform);
interpolator2->Initialize();
@ -960,16 +944,18 @@ void itkImageProcessor::GetProjectionImages(){
finalTransform->SetRotation(dtr * RZero[0], dtr * RZero[1], dtr * RZero[2]);
// The transform is determined by the parameters and the rotation center.
finalTransform->SetCenter(isocenter3D);
finalTransform->SetCenter(
m_3DProjectionOriginLPSZero);
using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
// The ResampleImageFilter is the driving force for the projection image generation.
ResampleFilterType::Pointer resampleFilter1 = ResampleFilterType::New();
resampleFilter1->SetInput(caster3D->GetOutput()); // Link the 3D volume.
resampleFilter1->SetInput(
image3DIn);
resampleFilter1->SetDefaultPixelValue(0);
resampleFilter1->SetNumberOfWorkUnits(1);
resampleFilter1->SetNumberOfWorkUnits(8);
// The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance
// have been set before registration. Here we only need to replace the initial
@ -997,7 +983,7 @@ void itkImageProcessor::GetProjectionImages(){
// Compute the origin (in mm) of the 2D Image
origin[0] = -image1res[0] * o2Dx;
origin[1] = -image1res[1] * o2Dy;
origin[2] = -scd;
origin[2] = -d_scd;
resampleFilter1->SetSize(size);
resampleFilter1->SetOutputSpacing(spacing);
@ -1009,11 +995,10 @@ void itkImageProcessor::GetProjectionImages(){
// Do the same thing for the output image 2.
ResampleFilterType::Pointer resampleFilter2 = ResampleFilterType::New();
resampleFilter2->SetInput(caster3D->GetOutput());
resampleFilter2->SetInput(
image3DIn);
resampleFilter2->SetDefaultPixelValue(0);
resampleFilter1->SetNumberOfWorkUnits(1);
resampleFilter2->SetNumberOfWorkUnits(8);
// The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance
@ -1037,7 +1022,7 @@ void itkImageProcessor::GetProjectionImages(){
// Compute the origin (in mm) of the 2D Image
origin[0] = -image2res[0] * o2Dx;
origin[1] = -image2res[1] * o2Dy;
origin[2] = -scd;
origin[2] = -d_scd;
resampleFilter2->SetSize(size);
resampleFilter2->SetOutputSpacing(spacing);
@ -1087,10 +1072,6 @@ void itkImageProcessor::GetProjectionImages(){
// InternalImageType::DirectionType LATtoHFS =
// this->CalcDRTImageDirections(projAngle2);
/* Again.. */
resampleFilter2->Update();

174
reg23Topograms/itkDTRrecon/itkImageProcessor.h
View File

@ -43,40 +43,40 @@ namespace itk
{
class CommandIterationUpdate : public itk::Command
{
public:
using Self = CommandIterationUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
//class CommandIterationUpdate : public itk::Command
//{
//public:
// using Self = CommandIterationUpdate;
// using Superclass = itk::Command;
// using Pointer = itk::SmartPointer<Self>;
// itkNewMacro(Self);
protected:
CommandIterationUpdate() = default;
//protected:
// CommandIterationUpdate() = default;
public:
using OptimizerType = itk::PowellOptimizer;
using OptimizerPointer = const OptimizerType *;
//public:
// using OptimizerType = itk::PowellOptimizer;
// using OptimizerPointer = const OptimizerType *;
void
Execute(itk::Object * caller, const itk::EventObject & event) override
{
Execute((const itk::Object *)caller, event);
}
// void
// Execute(itk::Object * caller, const itk::EventObject & event) override
// {
// Execute((const itk::Object *)caller, event);
// }
void
Execute(const itk::Object * object, const itk::EventObject & event) override
{
auto optimizer = dynamic_cast<OptimizerPointer>(object);
if (typeid(event) != typeid(itk::IterationEvent))
{
return;
}
std::cout << "Iteration: " << optimizer->GetCurrentIteration() << std::endl;
std::cout << "Similarity: " << optimizer->GetValue() << std::endl;
std::cout << "Position: " << optimizer->GetCurrentPosition() << std::endl;
}
};
// void
// Execute(const itk::Object * object, const itk::EventObject & event) override
// {
// auto optimizer = dynamic_cast<OptimizerPointer>(object);
// if (typeid(event) != typeid(itk::IterationEvent))
// {
// return;
// }
// std::cout << "Iteration: " << optimizer->GetCurrentIteration() << std::endl;
// std::cout << "Similarity: " << optimizer->GetValue() << std::endl;
// std::cout << "Position: " << optimizer->GetCurrentPosition() << std::endl;
// }
//};
/* reference string required for comparison with tag values */
@ -107,35 +107,42 @@ public:
/** Input data load methods*/
int load3DSerie(const char* );
int load2D(const char *);
int loadRTPlan(const char *);
int loadRTStruct(const char *);
/** Projection angles - Gantry angle */
/** Projection angles - Gantry angle IEC */
void SetProjectionAngle1(double);
void SetProjectionAngle2(double);
/** Source to axis distance - SAD */
/** Source to axis distance - SAD
* Often called source to object (SOD), (0018,1111) */
void SetSCD(double);
/** Intensity threshold used for image projection,
* any voxel below threshold is ignored */
* any voxel below threshold is ignored */
void SetIntensityThreshold(double);
/** Custom settings of the projection images */
void SetCustom_Isocenter(double,double,double);
void SetCustom_2Dres(double,double,double,double);
void SetCustom_2Dsize(int,int,int,int);
void SetCustom_2Dcenter(double dX1, double dY1, double dX2, double dY2);
/** Set transform parameters for 3D Volume */
void SetInitialTranslations(double, double, double);
void SetInitialRotations(double, double, double);
/** Get the Localizer intensity window and center for rendering */
double* GetImageIntensityWindow(int );
/** 2D-3D registration routines */
// int InitializeRegistration();
/** Initialize projection geometry */
void InitializeProjector();
// void RunRegistration();
/** Apply transform to CT volume.
* This moves the volume based on RT Plan info
* if those are available. */
void ApplyVolumeImportTransform();
/** Get the Localizer intensity window and
* center for rendering */
const double *GetImageIntensityWindow(int );
/** Compute Digital Localizers */
void GetProjectionImages();
@ -166,6 +173,7 @@ public:
typedef enum eImageOrientationType
{ NotDefined = 0, HFS = 1, FFS = 2}
tPatOrientation;
protected:
itkImageProcessor();
virtual ~itkImageProcessor();
@ -181,7 +189,7 @@ private:
void Finalize();
//void Finalize();
/** Image types */
@ -211,7 +219,6 @@ private:
using Input2DRescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, InternalImageType>;
using CastFilterType3D = itk::CastImageFilter<ImageType3D, InternalImageType>;
using ChangeInformationFilterType = itk::ChangeInformationImageFilter<InternalImageType>;
using ChangeInformationFilterInputType = itk::ChangeInformationImageFilter<ImageType3D>;
using ChangeInformationFilterInput2DType = itk::ChangeInformationImageFilter<ImageType2D>;
/** ITK to VTK filter */
@ -220,54 +227,44 @@ private:
TransformType::Pointer transform;
TransformType::Pointer
transform;
InterpolatorType::Pointer interpolator1;
InterpolatorType::Pointer interpolator2;
ImageReaderType3D::Pointer imageReader3D;
ImageSeriesReaderType3D::Pointer imageSeriesReader3D;
CommandIterationUpdate::Pointer observer;
CastFilterType3D::Pointer caster3D;
CastFilterType3D::Pointer caster2D1;
CastFilterType3D::Pointer caster2D2;
InterpolatorType::Pointer
interpolator1,
interpolator2;
TransformType::InputPointType isocenter3D;
ITKtoVTKFilterType::Pointer
toVTK2D1,
toVTK2D2,
toVTKLocalizer1,
toVTKLocalizer2;
ITKtoVTKFilterType::Pointer toVTK2D1;
ITKtoVTKFilterType::Pointer toVTK2D2;
ITKtoVTKFilterType::Pointer toVTKLocalizer1;
ITKtoVTKFilterType::Pointer toVTKLocalizer2;
InternalImageType::Pointer
image3DIn,
imageDRT1In,
imageDRT2In;
ImageType3D::Pointer image3DIn;
DuplicatorType::Pointer
m_LATSourceDupli,
m_PASourceDupli,
m_VolumeSourceDupli;
InternalImageType::Pointer imageDRT1In;
InternalImageType::Pointer imageDRT2In;
DuplicatorType::Pointer m_LATSourceDupli;
DuplicatorType::Pointer m_PASourceDupli;
ChangeInformationFilterInputType::Pointer
m_3DInputChangeInformation;
ChangeInformationFilterType::Pointer
m_3DInputChangeInformationToZero;
tPatOrientation
m_3DPatOrientation;
m_3DPatOrientation;
ImageType3D :: PointType
m_3DOriginLPS,
m_3DProjectionOriginLPS,
m_3DProjectionOriginLPSZero;
m_3DOriginLPS,
m_3DProjectionOriginLPS,
m_3DProjectionOriginLPSZero;
InternalImageType::DirectionType
LPStoIEC_Directions;
InternalImageType::DirectionType
CalcDRTImageDirections(double angle);
@ -276,6 +273,7 @@ private:
ImageType3D::PointType m_VolCOOV,
double dAngle
);
double
CalcProjectionAngleLPS(
tPatOrientation pOrient,
@ -288,33 +286,27 @@ private:
ImageType3D::PointType rtIsocenterLPS,
ImageType3D::PointType IEC2DCMMapT);
double image1res[2], image2res[2];
double image1center[2], image2center[2];
int image1Size[2], image2Size[2];
double isocenter[3];
double scd;
double threshold;
double
d_scd,
d_intensityThreshold;
double projAngle1, projAngle2;
double projAngle1_IEC, projAngle2_IEC;
double TZero[3], RZero[3];
bool customized_iso,customized_2DCX,customized_2DRES, customized_2DSIZE;
bool customized_iso,/*customized_2DCX,*/customized_2DRES, customized_2DSIZE;
bool image2D1Loaded, image2D2Loaded, image3DLoaded;
double image1IntensityWindow[2], image2IntensityWindow[2];
double
d_CTImportOffset_txiec,
d_CTImportOffset_tyiec,
d_CTImportOffset_tziec,
d_CTImportOffset_rxiec,
d_CTImportOffset_ryiec,
d_CTImportOffset_rziec;
d_image1IntensityWindow[2],
d_image2IntensityWindow[2];
InternalImageType::DirectionType Std_DRT2LPS;