458 lines
14 KiB
C++
458 lines
14 KiB
C++
|
|
/*
|
|
|
|
Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
|
|
version of incremental ray-tracing algorithm
|
|
|
|
gfattori 08.11.2021
|
|
|
|
*/
|
|
|
|
|
|
#ifndef ITKIMAGEPROCESSOR_H
|
|
#define ITKIMAGEPROCESSOR_H
|
|
|
|
|
|
#include "DRTMetaInformation.h"
|
|
|
|
#include "itkCommand.h"
|
|
|
|
#include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
|
|
#include "itkEuler3DTransform.h"
|
|
#include "itkResampleImageFilter.h"
|
|
#include "itkCastImageFilter.h"
|
|
#include "itkRescaleIntensityImageFilter.h"
|
|
#include "itkChangeInformationImageFilter.h"
|
|
#include "itkGDCMImageIO.h"
|
|
#include "itkMetaDataObject.h"
|
|
#include "itkImageDuplicator.h"
|
|
#include "itkObject.h"
|
|
#include "itkObjectFactory.h"
|
|
#include "itkSmartPointer.h"
|
|
|
|
#include "itkExhaustiveOptimizer.h"
|
|
#include "itkMutualInformationTwoImageToOneImageMetric.h"
|
|
#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h"
|
|
#include "itkPowellOptimizer.h"
|
|
#include "itkTwoProjectionImageRegistrationMethod.h"
|
|
#include "itkAmoebaOptimizer.h"
|
|
|
|
#include "itkImageToVTKImageFilter.h"
|
|
|
|
#include "itkImageFileReader.h"
|
|
#include "itkImageFileWriter.h"
|
|
#include "itkImageSeriesReader.h"
|
|
|
|
#include "gdcmGlobal.h"
|
|
|
|
#include "vtkImageData.h"
|
|
#include "vtkTransform.h"
|
|
#include "vtkSmartPointer.h"
|
|
|
|
#include "itkQtIterationUpdate.h"
|
|
|
|
namespace itk
|
|
{
|
|
|
|
class ITK_EXPORT itkImageProcessor : public itk::Object
|
|
{
|
|
|
|
constexpr static unsigned int Dimension = 3;
|
|
|
|
public:
|
|
/** Standard "Self" typedef. */
|
|
typedef itkImageProcessor Self;
|
|
/** Standard "Superclass" typedef. */
|
|
typedef Object Superclass;
|
|
/** Smart pointer typedef support */
|
|
typedef itk::SmartPointer<Self> Pointer;
|
|
typedef itk::SmartPointer<const Self> ConstPointer;
|
|
/** Method of creation through the object factory. */
|
|
itkNewMacro(Self);
|
|
/** Run-time type information (and related methods). */
|
|
itkTypeMacro(itkImageProcessor, Object);
|
|
|
|
|
|
CommandIterationUpdate::Pointer
|
|
optimizerObserver;
|
|
|
|
|
|
/** Input data load methods*/
|
|
int load3DSerieFromFolder(const char* );
|
|
int load3DSerieFromFiles( std::vector<std::string> );
|
|
int unload3DVolumeAndMeta();
|
|
|
|
int load2D(const char *);
|
|
int unload2DAndMeta(int);
|
|
int query2DimageReconstructionDiameter(const char*);
|
|
|
|
void loadRTPlanInfo(double, double, double, double, double ,double);
|
|
int unloadRTPlanAndMeta();
|
|
|
|
/** Projection angles - Gantry angle IEC */
|
|
void SetProjectionAngle1IEC(double);
|
|
void SetProjectionAngle2IEC(double);
|
|
void SetProjectionAngleOffsetIEC(double,double);
|
|
/** Get projection angles - Gantry angle LPS
|
|
* Only meaningful after a 3D Volume is loaded */
|
|
double GetProjectionAngle1LPS();
|
|
double GetProjectionAngle2LPS();
|
|
|
|
/** Source to axis distance - SAD
|
|
* Often called source to object (SOD), (0018,1111) */
|
|
void SetSCD(double,double);
|
|
void SetSCDOffset(double,double);
|
|
|
|
double GetSCD1();
|
|
double GetSCD2();
|
|
|
|
/** Panel offsets - panel offsets IEC */
|
|
void SetPanelOffsets(double,double);
|
|
/** Get projection angles - Gantry angle LPS
|
|
* Only meaningful after a 3D Volume is loaded */
|
|
double GetPanelOffset1();
|
|
double GetPanelOffset2();
|
|
|
|
/** Sets the degree of freedom for omptimzation*/
|
|
void SetDegreeOfFreedom(tDegreeOfFreedomEnum);
|
|
tDegreeOfFreedomEnum GetDegreeOfFreedom();
|
|
|
|
void SetCustom_ProjectionCenterOffsetFixedAxes_IEC(double,double,double,double,double,double);
|
|
|
|
void SetCustom_ProjectionCenterFixedAxes_IEC(double,double,double);
|
|
|
|
/** Intensity threshold used for image projection,
|
|
* any voxel below threshold is ignored */
|
|
void SetIntensityThreshold(double);
|
|
|
|
/** Custom settings of the projection images */
|
|
void SetCustom_2Dres(double,double,double,double);
|
|
void SetCustom_2Dsize(int,int,int,int);
|
|
void SetCustom_ImportTransform(double, double, double,
|
|
double, double, double);
|
|
void SetCustom_UpdateMetaInfo();
|
|
|
|
|
|
/** Set transform parameters for 3D Volume */
|
|
void SetInitialTranslations(double, double, double);
|
|
void SetInitialRotations(double, double, double);
|
|
|
|
/** Get transform parameters for 3D Volume */
|
|
void GetFinalTranslations(double&, double&, double&);
|
|
void GetFinalRotations(double&, double&, double&);
|
|
|
|
/** Set user transform parameters for 3D Volume in LPS*/
|
|
void SetUserTranslationsLPS(double, double, double);
|
|
void SetUserRotationsLPS(double, double, double);
|
|
|
|
/** Initialize the registration pipeline*/
|
|
void InitializeRegistration(double, double, eDegreeOfFreedomType);
|
|
/** Start the registration process*/
|
|
int StartRegistration(std::string extraInfo);
|
|
|
|
/** Get transform parameters for 3D Volume */
|
|
double* GetTransformParameters();
|
|
|
|
/** Initialize projection geometry */
|
|
void InitializeProjector();
|
|
|
|
/** Get Projection origin in LPS coordinates*/
|
|
const std::vector <double> GetLPStoProjectionGeoLPSOffset();
|
|
|
|
/** Get Projection origin in LPS coordinates*/
|
|
const std::vector <double> GetRTImportOffset();
|
|
|
|
/** Get the Localizer intensity window and
|
|
* center for rendering */
|
|
double GetLocalizerDisplayWindowLevel(int );
|
|
double GetLocalizerDisplayWindowWidth(int );
|
|
|
|
/** Compute Digital Localizers */
|
|
void GetProjectionImages();
|
|
|
|
/** Conveniency methods to get VTK images for rendering */
|
|
vtkImageData* GetProjection1VTK();
|
|
vtkImageData* GetProjection2VTK();
|
|
|
|
vtkImageData* GetLocalizer1VTK();
|
|
vtkImageData* GetLocalizer2VTK();
|
|
|
|
vtkMatrix4x4 * GetProjection1VTKTransform();
|
|
vtkMatrix4x4 * GetProjection2VTKTransform();
|
|
|
|
/** Debug writers */
|
|
void WriteProjectionImages();
|
|
void Write2DImages();
|
|
|
|
/** Get the current cost function value from the optimizer*/
|
|
double GetOptimizerValue();
|
|
|
|
/** Get the cost function value for the current transform*/
|
|
double GetCurrentMetricValue();
|
|
|
|
void SetROI(double, double, double, double);
|
|
|
|
void SetOptimizer(std::string);
|
|
void SetMetric(std::string);
|
|
void SetFullROI(bool);
|
|
void SetMaxNumberOfIterations(int);
|
|
void SetRegionFixed1(int,int,int,int);
|
|
void SetRegionFixed2(int,int,int,int);
|
|
|
|
/** Optimizer which tries to find the minimn (Powell Optimizer)*/
|
|
using OptimizerType = itk::PowellOptimizer;
|
|
|
|
protected:
|
|
/** Various pixel types */
|
|
using InternalPixelType = float;
|
|
using PixelType2D = short;
|
|
using PixelType3D = short;
|
|
using OutputPixelType = unsigned short;
|
|
|
|
using ImageType3D = itk::Image<PixelType3D, Dimension>;
|
|
using InternalImageType = itk::Image<InternalPixelType, Dimension>;
|
|
|
|
itkImageProcessor();
|
|
virtual ~itkImageProcessor();
|
|
void PrintSelf(std::ostream& os, itk::Indent indent) const;
|
|
|
|
private:
|
|
itkImageProcessor(const Self&); //purposely not implemented
|
|
void operator=(const Self&); //purposely not implemented
|
|
|
|
/** Fill Meta after 3D volume load */
|
|
int fill3DVolumeMeta(InternalImageType::Pointer,
|
|
tPatOrientation);
|
|
|
|
/** Image types */
|
|
using ImageType2D = itk::Image<PixelType3D, Dimension>;
|
|
|
|
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
|
|
|
|
/** The following lines define each of the components used in the
|
|
registration: The transform, optimizer, metric, interpolator and
|
|
the registration method itself. */
|
|
using TransformType = itk::Euler3DTransform<double>;
|
|
using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
|
|
using ResampleFilterType = itk::ResampleImageFilter<InternalImageType, InternalImageType>;
|
|
|
|
/** Optimizer which tries to find the minimn (Powell Optimizer)*/
|
|
using AmoebaOptimizerType = itk::AmoebaOptimizer;
|
|
/** Optimizer which samples the whol space */
|
|
using ExhaustiveOptimizerType = itk::ExhaustiveOptimizer;
|
|
/** Metric to calculate similarites between images*/
|
|
using MetricType = itk::NormalizedCorrelationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
|
|
using MIMetricType = itk::MutualInformationTwoImageToOneImageMetric<InternalImageType, InternalImageType>;
|
|
/** The thing which actuall does the image registration*/
|
|
using RegistrationType = itk::TwoProjectionImageRegistrationMethod<InternalImageType, InternalImageType>;
|
|
|
|
using RoiForRegistration = itk::ImageRegion<Dimension>;
|
|
|
|
/** Image reader types */
|
|
using ImageReaderType2D = itk::ImageFileReader<ImageType2D>;
|
|
using ImageReaderType3D = itk::ImageFileReader<ImageType3D>;
|
|
using ImageSeriesReaderType3D = itk::ImageSeriesReader<ImageType3D>;
|
|
using ImageIOType = itk::GDCMImageIO;
|
|
using DuplicatorType = itk::ImageDuplicator<InternalImageType>;
|
|
|
|
/** widely used filters: flip, cast, rescale, ... */
|
|
using Input2DRescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, InternalImageType>;
|
|
using CastFilterType3D = itk::CastImageFilter<ImageType3D, InternalImageType>;
|
|
using ChangeInformationFilterType = itk::ChangeInformationImageFilter<InternalImageType>;
|
|
using ChangeInformationFilterInput2DType = itk::ChangeInformationImageFilter<ImageType2D>;
|
|
|
|
/** ITK to VTK filter */
|
|
using ITKtoVTKFilterType = itk::ImageToVTKImageFilter<OutputImageType>;
|
|
|
|
void
|
|
CalculateExternalUserTransform(
|
|
TransformType::Pointer transform,
|
|
DRTImageMetaInformation::Pointer imageMetaInfo);
|
|
|
|
// TransformType::Pointer
|
|
// CalculateInternalTransform(
|
|
// ImageType3D::PointType m_TranslationOffset,
|
|
// ImageType3D::PointType m_RotationOffset,
|
|
// ImageType3D::PointType m_TranslationUser,
|
|
// ImageType3D::PointType m_RotationUser,
|
|
// ImageType3D::PointType m_ProjectionTransformCenter,
|
|
// ImageType3D::PointType m_UserTransformCenter,
|
|
// ImageType3D::PointType m_OffsetTransformCenter,
|
|
// InternalImageType::DirectionType m_IECtoLPSDirections
|
|
// );
|
|
|
|
/* Calculate the transform used in siddon.
|
|
* The isocentric transform is mapped to the calibrated center of projection */
|
|
TransformType::Pointer
|
|
CalculateInternalTransformV2(
|
|
ImageType3D::PointType m_TranslationOffset,
|
|
ImageType3D::PointType m_RotationOffset,
|
|
ImageType3D::PointType m_TranslationUser,
|
|
ImageType3D::PointType m_RotationUser,
|
|
ImageType3D::PointType m_CalibratedProjectionCenter,
|
|
ImageType3D::PointType m_RTIsocenter,
|
|
InternalImageType::DirectionType m_IECtoLPSDirections
|
|
);
|
|
|
|
|
|
|
|
TransformType::Pointer
|
|
transform1,
|
|
transform2;
|
|
|
|
InterpolatorType::Pointer
|
|
interpolator1,
|
|
interpolator2;
|
|
|
|
ITKtoVTKFilterType::Pointer
|
|
toVTK2D1,
|
|
toVTK2D2,
|
|
toVTKLocalizer1,
|
|
toVTKLocalizer2;
|
|
|
|
InternalImageType::Pointer
|
|
image3DIn,
|
|
imageDRT1In,
|
|
imageDRT2In;
|
|
|
|
RegistrationType::Pointer
|
|
registration;
|
|
MetricType::Pointer
|
|
metric;
|
|
MIMetricType::Pointer
|
|
mimetric;
|
|
OptimizerType::Pointer
|
|
optimizer;
|
|
AmoebaOptimizerType::Pointer
|
|
amoebaoptimizer;
|
|
ExhaustiveOptimizerType::Pointer
|
|
exhaustiveOptimizer;
|
|
ExhaustiveCommandIterationUpdate::Pointer
|
|
exhaustiveOptimizerObserver;
|
|
|
|
DuplicatorType::Pointer
|
|
m_LATSourceDupli,
|
|
m_PASourceDupli,
|
|
m_VolumeSourceDupli;
|
|
|
|
ChangeInformationFilterType::Pointer
|
|
m_3DInputChangeInformationToZero;
|
|
|
|
/** Apply transform to CT volume.
|
|
* This moves the volume based on RT Plan info
|
|
* if those are available. */
|
|
void UpdateProjectionGeometryMeta();
|
|
|
|
/* The following functions do not rely on class variables,
|
|
* Only input variables are used... */
|
|
|
|
InternalImageType::DirectionType
|
|
CalcDRTImageDirections(double angle,
|
|
InternalImageType::DirectionType DRT2LPS);
|
|
|
|
ImageType3D::PointType
|
|
CalcDRTImageOrigin(
|
|
ImageType3D::PointType m_DRTOrigin,
|
|
ImageType3D::PointType m_DRTCOV,
|
|
ImageType3D::PointType m_ProjectionCenterLPS,
|
|
double dAngle,
|
|
InternalImageType::DirectionType stdDRT2LPS
|
|
);
|
|
|
|
|
|
// ImageType3D::PointType
|
|
// CalcDRTImageOffset(
|
|
// ImageType3D::PointType m_DRTOffset,
|
|
// double dAngle,
|
|
// InternalImageType::DirectionType stdDRT2LPS
|
|
// );
|
|
|
|
TransformType::Pointer
|
|
MapTransformToNewOrigin(
|
|
ImageType3D::PointType m_COR,
|
|
ImageType3D::PointType m_Translations,
|
|
ImageType3D::PointType m_Rotations
|
|
);
|
|
|
|
double
|
|
CalcProjectionAngleLPS(
|
|
tPatOrientation pOrient,
|
|
double pAngleIEC);
|
|
|
|
/** Apply transform to CT volume.
|
|
* This moves the volume based on RT Plan info
|
|
* if those are available. */
|
|
ImageType3D::PointType
|
|
CalcImportVolumeOffset(
|
|
ImageType3D::PointType rtCouchOffset,
|
|
InternalImageType::DirectionType VolumeLPStoIEC,
|
|
ImageType3D::PointType rtIsocenterLPS,
|
|
ImageType3D::PointType IEC2DCMMapT,
|
|
ImageType3D::PointType IEC2DCMMapR);
|
|
|
|
|
|
// The ResampleImageFilter is the driving force for the projection image generation.
|
|
ResampleFilterType::Pointer
|
|
resampleFilter1,
|
|
resampleFilter2;
|
|
|
|
ChangeInformationFilterType::Pointer
|
|
filter1,
|
|
filter2;
|
|
|
|
InternalImageType::DirectionType
|
|
Std_DRT2LPS;
|
|
|
|
vtkMatrix4x4
|
|
* m_Projection1VTKTransform,
|
|
* m_Projection2VTKTransform;
|
|
|
|
/*Transformation Parameters */
|
|
double dTransfParam[6];
|
|
|
|
/**
|
|
* Many meta containers
|
|
*/
|
|
CTVolumeImageMetaInformation::Pointer
|
|
m_CTMetaInfo;
|
|
|
|
DRTProjectionGeometryImageMetaInformation::Pointer
|
|
m_DRTGeometryMetaInfo;
|
|
|
|
DRTImageMetaInformation::Pointer
|
|
m_DRTImage1MetaInfo,
|
|
m_DRTImage2MetaInfo;
|
|
|
|
TopogramImageMetaInformation::Pointer
|
|
m_TImage1MetaInfo,
|
|
m_TImage2MetaInfo;
|
|
|
|
RTGeometryMetaInformation::Pointer
|
|
m_RTMetaInfo;
|
|
|
|
TransformMetaInformation::Pointer
|
|
m_TransformMetaInfo;
|
|
|
|
|
|
double m_OptmizerValue;
|
|
int m_MaxNumberOfIterations;
|
|
|
|
RoiForRegistration
|
|
roiAutoReg1,
|
|
roiAutoReg2;
|
|
|
|
|
|
bool m_UseExhaustiveOptmizer;
|
|
bool m_UseAmeobaOptimizer;
|
|
bool m_UseFullROI;
|
|
bool m_UseMutualInformation;
|
|
bool m_UseDumptoFile;
|
|
|
|
};
|
|
|
|
|
|
}
|
|
|
|
|
|
#endif
|