fattori_g/reg23Topograms
0
0
Fork 0
Code Issues 1 Pull Requests Actions Packages Projects Releases Wiki Activity
Files
master
reg23Topograms/itkReg23DRT/itkImageProcessor.cpp
Giovanni Fattori 2361354982 A bit of code comments cleanup
2025-05-14 23:20:02 +02:00

2408 lines
69 KiB
C++
Raw Permalink Blame History

/*
Calculate Ditigally Reconstructed Topogram from a CT dataset using a modified
version of incremental ray-tracing algorithm
gfattori 08.11.2021
*/
#include "itkImageProcessor.h"
#include "itkCreateObjectFunction.h"
#include "itkDRTHelpers.h"
#include "itkVersion.h"
#include <gdcmReader.h>
#include <itkMatrix.h>
#include "itkImageDuplicator.h"
#include "itkPermuteAxesImageFilter.h"
#include "itkMatrix.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkOrientImageFilter.h"
#include <gdcmIPPSorter.h>
#include "itkCommand.h"
#include "itkTimeProbesCollectorBase.h"
#include <string>
namespace itk
{
itkImageProcessor::itkImageProcessor()
{
iNWUnits = 1;
transform1 = TransformType::New();
transform1->SetIdentity();
transform2 = TransformType::New();
transform2->SetIdentity();
m_InternalTransf1 = InternalTransformMetaInformation::New();
m_InternalTransf2 = InternalTransformMetaInformation::New();
interpolator1 = InterpolatorType::New();
interpolator2 = InterpolatorType::New();
toVTK2D1 = ITKtoVTKFilterType::New();
toVTK2D2 = ITKtoVTKFilterType::New();
toVTKLocalizer1 = ITKtoVTKFilterType::New();
toVTKLocalizer2 = ITKtoVTKFilterType::New();
m_LATSourceDupli = DuplicatorType::New();
m_PASourceDupli = DuplicatorType::New();
m_VolumeSourceDupli = DuplicatorType::New();
Std_DRT2LPS.SetIdentity();
for(int iIdx = 0 ; iIdx < 3; iIdx++){
for(int jIdx = 0 ; jIdx < 3; jIdx++){
Std_DRT2LPS.GetVnlMatrix()[iIdx][jIdx] = Standard_DRT2LPS [jIdx+iIdx*3];
}
}
m_Projection1VTKTransform = vtkMatrix4x4::New();
m_Projection1VTKTransform->Identity();
m_Projection2VTKTransform = vtkMatrix4x4::New();
m_Projection2VTKTransform->Identity();
/* Set to NULL the metainfo */
m_CTMetaInfo = NULL;
m_TImage1MetaInfo = NULL;
m_TImage2MetaInfo = NULL;
m_DRTImage1MetaInfo = NULL;
m_DRTImage2MetaInfo = NULL;
m_RTMetaInfo = NULL;
m_TransformMetaInfo = NULL;
m_TransformMetaInfo = TransformMetaInformation::New();
m_r23MetaInfo = NULL;
m_r23MetaInfo = R23MetaInformation::New();
m_PowellOptimizerMetaInfo = NULL;
m_PowellOptimizerMetaInfo = PowellOptimizerMetaInformation::New();
m_AmoebaOptimizerMetaInfo = NULL;
m_AmoebaOptimizerMetaInfo = AmoebaOptimizerMetaInformation::New();
m_MIMetricMetaInfo = NULL;
m_MIMetricMetaInfo = MIMetricMetaInformation::New();
m_NCCMetricMetaInfo = NULL;
m_NCCMetricMetaInfo = NCCMetricMetaInformation::New();
/* Initialise the projection geoemtry with defaults */
m_DRTGeometryMetaInfo = NULL;
m_DRTGeometryMetaInfo
= DRTProjectionGeometryImageMetaInformation::New();
m_DRTGeometryMetaInfo->SetSCD1(570.);
m_DRTGeometryMetaInfo->SetSCD2(570.);
m_DRTGeometryMetaInfo->SetProjectionAngle1IEC(180.);
m_DRTGeometryMetaInfo->SetProjectionAngle2IEC(90.);
m_DRTGeometryMetaInfo->SetIntensityThreshold(-300.);
ImageType3D::PointType
Point3D(3);
Point3D[0]=0.;
Point3D[1]=0.;
Point3D[2]=175.;
m_DRTGeometryMetaInfo->SetProjectionCenter(Point3D);
ImageType3D::SizeType
ImageSize;
ImageSize[0]=512;
ImageSize[1]=512;
ImageSize[2]=1;
m_DRTGeometryMetaInfo->SetDRT1Size(ImageSize);
m_DRTGeometryMetaInfo->SetDRT2Size(ImageSize);
ImageType3D::SpacingType
ImageRes(3);
ImageRes[0] = 1.;
ImageRes[1] = 1.;
ImageRes[2] = 1.;
m_DRTGeometryMetaInfo->SetDRT1Spacing(ImageRes);
m_DRTGeometryMetaInfo->SetDRT2Spacing(ImageRes);
Point3D[0]=0.;
Point3D[1]=0.;
Point3D[2]=0.;
m_DRTGeometryMetaInfo->SetIECS2IECScannerT(Point3D);
m_DRTGeometryMetaInfo->SetIECS2IECScannerR(Point3D);
m_DRTGeometryMetaInfo->SetProjectionCenterOffset1(Point3D);
m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(Point3D);
m_DRTGeometryMetaInfo->SetDegreeOfFreedom(eDegreeOfFreedomType::THREE_DOF);
m_DRTGeometryMetaInfo->SetHandleRotationImpOffset(eHandlingRotImpTransform::ALWAYS_USE);
optimizerObserver = CommandIterationUpdate::New();
ExhaustiveOptimizerObserver = ExhaustiveCommandIterationUpdate::New();
m_R23 = itkReg23::New();
m_R23->SetOptimizerObserver(optimizerObserver);
}
itkImageProcessor::~itkImageProcessor()
{
}
void itkImageProcessor::SetNumberOfWorkingUnits(int iN){
if(iN>0 && iN<50){
iNWUnits = iN;
}
}
const CTVolumeImageMetaInformation::Pointer
itkImageProcessor::GetCTMetaInfo(
){
return m_CTMetaInfo;
}
const RTGeometryMetaInformation::Pointer
itkImageProcessor::GetRTMetaInfo(
){
return m_RTMetaInfo;
}
const DRTProjectionGeometryImageMetaInformation::Pointer
itkImageProcessor::GetDRTGeoMetaInfo(){
return m_DRTGeometryMetaInfo;
}
void itkImageProcessor::PrintSelf(std::ostream& os, Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
void itkImageProcessor::SetIntensityThreshold(double dT)
{
m_DRTGeometryMetaInfo->SetIntensityThreshold(dT);
}
void itkImageProcessor::SetSCDOffset(double dOff1, double dOff2)
{
m_DRTGeometryMetaInfo->SetSCD1Offset(dOff1);
m_DRTGeometryMetaInfo->SetSCD2Offset(dOff2);
}
void itkImageProcessor::SetSCD(double dDist1, double dDist2)
{
m_DRTGeometryMetaInfo->SetSCD1(dDist1);
m_DRTGeometryMetaInfo->SetSCD2(dDist2);
}
double itkImageProcessor::GetSCD1(){
return
m_DRTImage1MetaInfo->GetSCD();
}
double itkImageProcessor::GetSCD2(){
return
m_DRTImage2MetaInfo ->GetSCD();
}
double itkImageProcessor::GetPanelOffsetPGeo(tProjOrientationType ImgPrj){
if(m_CTMetaInfo == NULL ||
m_DRTGeometryMetaInfo == NULL)
return 0.;
switch (ImgPrj) {
case tProjOrientationType::PA:
return
this->CalcPanelOffsetPGeo(m_CTMetaInfo->GetPatientOrientation(),
-m_DRTGeometryMetaInfo->GetPanel1Offset());
break;
case tProjOrientationType::LAT:
return
this->CalcPanelOffsetPGeo(m_CTMetaInfo->GetPatientOrientation(),
-m_DRTGeometryMetaInfo->GetPanel2Offset());
break;
case tProjOrientationType::NA:
return 0.;
break;
default:
return 0.;
break;
}
}
void itkImageProcessor::SetPanelOffsets(double dOff1, double dOff2)
{
m_DRTGeometryMetaInfo->SetPanel1Offset(dOff1);
m_DRTGeometryMetaInfo->SetPanel2Offset(dOff2);
}
double itkImageProcessor::GetPanelOffset1(){
return
m_DRTGeometryMetaInfo->GetPanel1Offset();
}
double itkImageProcessor::GetPanelOffset2(){
return
m_DRTGeometryMetaInfo ->GetPanel2Offset();
}
tDegreeOfFreedomEnum
itkImageProcessor::GetDegreeOfFreedom()
{
return
m_r23MetaInfo->GetDegreeOfFreedom();
}
void itkImageProcessor::SetUseRotationsForImportOffset(bool bVal){
if(m_DRTGeometryMetaInfo == NULL){
return;
}
m_DRTGeometryMetaInfo->SetUseRotationsForHiddenTransform(bVal);
return;
}
void itkImageProcessor::SetCustom_ImportTransform(double dTx,
double dTy,
double dTz,
double dRx,
double dRy,
double dRz)
{
if(m_DRTGeometryMetaInfo == NULL){
//todo
}
ImageType3D::PointType
Punto;
Punto[0] = dTx;
Punto[1] = dTy;
Punto[2] = dTz;
m_DRTGeometryMetaInfo->SetIECS2IECScannerT(Punto);
Punto[0] = dRx;
Punto[1] = dRy;
Punto[2] = dRz;
m_DRTGeometryMetaInfo->SetIECS2IECScannerR(Punto);
}
void itkImageProcessor::SetCustom_ProjectionCenterOffsetFixedAxes_IEC(
double dX1,
double dY1,
double dZ1,
double dX2,
double dY2,
double dZ2)
{
typedef itk::Point<double, 3> PointType;
PointType m_point;
m_point [0] = dX1;
m_point [1] = dY1;
m_point [2] = dZ1;
m_DRTGeometryMetaInfo->SetProjectionCenterOffset1(m_point);
m_point [0] = dX2;
m_point [1] = dY2;
m_point [2] = dZ2;
m_DRTGeometryMetaInfo->SetProjectionCenterOffset2(m_point);
}
void itkImageProcessor::SetCustom_ProjectionCenterFixedAxes_IEC(
double dX1,
double dY1,
double dZ1)
{
typedef itk::Point<double, 3> PointType;
PointType m_point;
m_point [0] = dX1;
m_point [1] = dY1;
m_point [2] = dZ1;
m_DRTGeometryMetaInfo->SetProjectionCenter(m_point);
}
void itkImageProcessor::SetCustom_2Dres(double nX1,double nY1,double nX2,double nY2)
{
if(m_DRTGeometryMetaInfo == NULL) {
// todo
}
ImageType3D::SpacingType Spacing;
Spacing [0] = nX1;
Spacing [1] = nY1;
Spacing [2] = 1.;
m_DRTGeometryMetaInfo->SetDRT1Spacing(Spacing);
Spacing [0] = nX2;
Spacing [1] = nY2;
Spacing [2] = 1.;
m_DRTGeometryMetaInfo->SetDRT2Spacing(Spacing);
}
void itkImageProcessor::SetCustom_2Dsize(int nX1, int nY1,int nX2,int nY2)
{
if(m_DRTGeometryMetaInfo == NULL) {
// todo
}
ImageType3D::SizeType Size;
Size [0] = nX1;
Size [1] = nY1;
Size [2] = 1.;
m_DRTGeometryMetaInfo->SetDRT1Size(Size);
Size [0] = nX2;
Size [1] = nY2;
Size [2] = 1.;
m_DRTGeometryMetaInfo->SetDRT2Size(Size);
}
void itkImageProcessor::SetCustom_UpdateMetaInfo(){
this->UpdateProjectionGeometryMeta();
}
int itkImageProcessor::unload3DVolumeAndMeta(){
std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
if (m_VolumeSourceDupli) {
m_VolumeSourceDupli = NULL;
m_VolumeSourceDupli = DuplicatorType::New();
}
m_CTMetaInfo = NULL;
m_DRTImage1MetaInfo = NULL;
m_DRTImage2MetaInfo = NULL;
m_TransformMetaInfo = NULL;
m_TransformMetaInfo = TransformMetaInformation::New();
this->ResetROIRegions();
return 1;
}
int itkImageProcessor::load3DSerieFromFiles( std::vector<std::string> v_fnames){
tPatOrientation m_PatOrientation
= tPatOrientation::NotDefined;
/* Since are not sure what we get here,
* we run the IPP sorter. */
using IppSorterType = gdcm::IPPSorter;
IppSorterType ZSorter;
ZSorter.SetComputeZSpacing( true );
ZSorter.SetZSpacingTolerance( 1e-3 );
bool b = ZSorter.Sort(v_fnames);
if( !b )
{
std::cerr << "itkImageProcessor::load3DSerieFromFiles :: Failed to sort files based on IPP." << std::endl;
return EXIT_FAILURE;
}
const std::vector<std::string> & v_sortedFnames =
ZSorter.GetFilenames();
try
{
using ImageIOType = itk::GDCMImageIO;
ImageIOType::Pointer dicomIO = ImageIOType::New();
ImageSeriesReaderType3D::Pointer
imageSeriesReader3D = ImageSeriesReaderType3D::New();
imageSeriesReader3D->SetImageIO(dicomIO);
imageSeriesReader3D->SetFileNames(v_sortedFnames);
imageSeriesReader3D->SetNumberOfWorkUnits(iNWUnits);
imageSeriesReader3D->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt
imageSeriesReader3D->Update();
if(v_fnames.size()>0){
gdcm::Reader R;
R.SetFileName(v_sortedFnames.at(0).c_str());
if(!R.Read())
{
cerr<<"ERROR: cannot read file: "<<v_sortedFnames.at(0).c_str()<<endl;
return -1;
}
const gdcm::File &file = R.GetFile();
const gdcm::DataSet &ds = file.GetDataSet();
char* sTmpString = new char [255];
strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0018,0x5100), ds));
*std::remove(sTmpString, sTmpString + strlen(sTmpString), ' ') = 0;
if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::HFS])){
m_PatOrientation = eImageOrientationType::HFS;
} else if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::FFS])){
m_PatOrientation = eImageOrientationType::FFS;
} else {
m_PatOrientation = eImageOrientationType::NotDefined;
}
delete [] (sTmpString);
} else {
return -1;
}
CastFilterType3D::Pointer caster3D =
CastFilterType3D::New();
if(m_ResamplepCT == false){
caster3D->SetInput(imageSeriesReader3D->GetOutput());
caster3D->Update();
}else{
std::cout << "Resampling pCT @ spacing in Z " << m_SamplingLNG << std::endl;
LinearInterpolatorType::Pointer
linearInterpolator = LinearInterpolatorType::New();
// Set identity transform
TransformType::Pointer
transform = TransformType::New();
transform->SetIdentity();
const typename ImageType3D::SpacingType& inputSpacing =
imageSeriesReader3D->GetOutput()->GetSpacing();
const typename ImageType3D::RegionType& inputRegion =
imageSeriesReader3D->GetOutput()->GetLargestPossibleRegion();
const typename ImageType3D::SizeType& inputSize =
inputRegion.GetSize();
ImageType3D::SpacingType outputSpacing;
outputSpacing[0] = inputSpacing[0];
outputSpacing[1] = inputSpacing[1];
outputSpacing[2] = m_SamplingLNG;
ImageType3D::SizeType outputSize;
outputSize[0] = static_cast<SizeValueType>(
inputSize[0] * inputSpacing[0] / outputSpacing[0] + .5);
outputSize[1] = static_cast<SizeValueType>(
inputSize[1] * inputSpacing[1] / outputSpacing[1] + .5);
outputSize[2] = static_cast<SizeValueType>(
inputSize[2] * inputSpacing[2] / outputSpacing[2] + .5);
ResampleInputFilterType::Pointer
resampler = ResampleInputFilterType::New();
resampler->SetInput( imageSeriesReader3D->GetOutput() );
resampler->SetTransform( transform );
resampler->SetInterpolator( linearInterpolator );
resampler->SetDefaultPixelValue(-1024);
resampler->SetOutputOrigin( imageSeriesReader3D->GetOutput()->GetOrigin() );
resampler->SetOutputSpacing( outputSpacing );
resampler->SetOutputDirection( imageSeriesReader3D->GetOutput()->GetDirection() );
resampler->SetSize( outputSize );
resampler->Update();
caster3D->SetInput(resampler->GetOutput());
caster3D->Update();
}
if (m_VolumeSourceDupli == NULL)
std::cout << "NEVER HERE m_VolumeSourceDupli" << std::endl;
m_VolumeSourceDupli->SetInputImage(caster3D->GetOutput());
m_VolumeSourceDupli->Update();
caster3D = NULL;
}
catch (itk::ExceptionObject & ex)
{
// std::cout << ex << std::endl;
return EXIT_FAILURE;
}
InternalImageType::Pointer m_InputImage =
m_VolumeSourceDupli->GetOutput();
return
this->fill3DVolumeMeta(m_InputImage, m_PatOrientation);
}
/** Fill Meta after 3D volume load */
int itkImageProcessor::fill3DVolumeMeta(
InternalImageType::Pointer m_InputImage,
tPatOrientation m_PatOrientation){
if(m_CTMetaInfo != NULL){
std::cout << "NEVER HERE m_CTMetaInfo" << std::endl;
m_CTMetaInfo = NULL;
}
if(m_DRTImage1MetaInfo != NULL){
std::cout << "NEVER HERE m_DRTImage1MetaInfo" << std::endl;
m_DRTImage1MetaInfo = NULL;
}
if(m_DRTImage2MetaInfo != NULL){
std::cout << "NEVER HERE m_DRTImage2MetaInfo" << std::endl;
m_DRTImage2MetaInfo = NULL;
}
if(m_RTMetaInfo != NULL){
std::cout << "NEVER HERE m_RTMetaInfo" << std::endl;
m_RTMetaInfo = NULL;
}
/* copy useful meta information into the CT container */
m_CTMetaInfo = CTVolumeImageMetaInformation::New();
m_CTMetaInfo->SetPatientOrientation(m_PatOrientation);
m_CTMetaInfo->SetSpacing(m_InputImage->GetSpacing());
m_CTMetaInfo->SetSize(
m_InputImage->GetBufferedRegion().GetSize() );
std::cout<<"ORIGIN: "<<m_InputImage->GetOrigin()<<std::endl;
m_CTMetaInfo->SetOriginLPS(m_InputImage->GetOrigin());
m_CTMetaInfo->SetImageDirections(m_InputImage->GetDirection());
ImageType3D::PointType
PointOffset;
PointOffset.Fill(0.);
m_CTMetaInfo->SetImportOffset(PointOffset);
/* initialise DRT meta */
m_DRTImage1MetaInfo = DRTImageMetaInformation::New();
m_DRTImage1MetaInfo->SetProjectionAngleLPS(
this->CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() +
m_DRTGeometryMetaInfo->GetProjectionAngle1OffsetIEC())
);
m_DRTImage1MetaInfo->SetSCD(
m_DRTGeometryMetaInfo->GetSCD1() +
m_DRTGeometryMetaInfo->GetSCD1Offset()
);
/* Calculate projection angle IEC to LPS */
m_DRTImage2MetaInfo = DRTImageMetaInformation::New();
m_DRTImage2MetaInfo->SetProjectionAngleLPS(
this->CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle2IEC() +
m_DRTGeometryMetaInfo->GetProjectionAngle2OffsetIEC())
);
m_DRTImage2MetaInfo->SetSCD(
m_DRTGeometryMetaInfo->GetSCD2()+
m_DRTGeometryMetaInfo->GetSCD2Offset()
);
std::cout<< "ImportOffset"<< m_CTMetaInfo->GetImportOffset() <<std::endl;
this->UpdateProjectionGeometryMeta();
return EXIT_SUCCESS;
}
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
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
itkImageProcessor::CalcPanelOffsetPGeo(
tPatOrientation pOrient,
double dPOffset){
double dVal = dPOffset;//0.;
switch (pOrient) {
case tPatOrientation :: HFS:
break;
case tPatOrientation :: FFS:
dVal *= -1;
break;
default:
break;
}
return
dVal;
}
double
itkImageProcessor::CalcProjectionAngleLPS(
tPatOrientation pOrient,
double pAngleIEC){
double dProj = pAngleIEC;
InternalImageType::DirectionType
currIECtoLPS;
/* NOTE that we transpose on the fly... */
switch (pOrient) {
case tPatOrientation :: HFS:
for(int iIdx = 0 ; iIdx < 3; iIdx++){
for(int jIdx = 0 ; jIdx < 3; jIdx++){
currIECtoLPS.GetVnlMatrix()[jIdx][iIdx] = HFS2IEC [jIdx+iIdx*3];
}
}
break;
case tPatOrientation :: FFS:
for(int iIdx = 0 ; iIdx < 3; iIdx++){
for(int jIdx = 0 ; jIdx < 3; jIdx++){
currIECtoLPS.GetVnlMatrix()[jIdx][iIdx] = FFS2IEC [jIdx+iIdx*3];
}
}
break;
default:
break;
}
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);
InternalImageType::DirectionType
DRTtoCurrPatientOrient;
DRTtoCurrPatientOrient =
currIECtoLPS * DRTDirsIEC;
ImageType3D::PointType SourcePositionIEC;
SourcePositionIEC [0] = 0;
SourcePositionIEC [1] = 0;
SourcePositionIEC [2] = 1;
ImageType3D::PointType SourcePositionPatOrient =
DRTtoCurrPatientOrient * SourcePositionIEC;
dProj =
(atan2(SourcePositionPatOrient [1], SourcePositionPatOrient [0]) - atan2(-1,0)) * 180 / itk::Math::pi;
return
dProj;
}
void itkImageProcessor::SetProjectionAngleOffsetIEC(double dOff1, double dOff2)
{
m_DRTGeometryMetaInfo->SetProjectionAngle1OffsetIEC(dOff1);
m_DRTGeometryMetaInfo->SetProjectionAngle2OffsetIEC(dOff2);
}
void itkImageProcessor::SetProjectionAngle1IEC(double dGantryAngle)
{
m_DRTGeometryMetaInfo->SetProjectionAngle1IEC(dGantryAngle);
}
double itkImageProcessor::GetProjectionAngle1LPS(){
return
m_DRTImage1MetaInfo->GetProjectionAngleLPS();
}
void itkImageProcessor::SetProjectionAngle2IEC(double dGantryAngle)
{
m_DRTGeometryMetaInfo->SetProjectionAngle2IEC(dGantryAngle);
}
double itkImageProcessor::GetProjectionAngle2LPS(){
return
m_DRTImage2MetaInfo->GetProjectionAngleLPS();
}
int itkImageProcessor::unload2DAndMeta(int iImgType){
switch (iImgType) {
case eProjectionOrientationType::PA:
std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << " PA "<< std::endl;
m_TImage1MetaInfo = NULL;
if (m_PASourceDupli) {
m_PASourceDupli = NULL;
m_PASourceDupli = DuplicatorType::New();
}
break;
case eProjectionOrientationType::LAT:
std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << " LAT "<< std::endl;
m_TImage2MetaInfo = NULL;
if (m_LATSourceDupli) {
m_LATSourceDupli = NULL;
m_LATSourceDupli = DuplicatorType::New();
}
break;
}
return 1;
}
int itkImageProcessor::load2D(const char * pcFName){
/* Check if we can open the file */
gdcm::Reader R;
R.SetFileName(pcFName);
if(!R.Read())
{ cerr<<"ERROR: cannot read file: "<<pcFName<<endl;
return -1;
}
/* Check if it's a localizer image */
const gdcm::File &file = R.GetFile();
const gdcm::DataSet &ds = file.GetDataSet();
char* sTmpString = new char [255];
strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0008,0x0008), ds));
std::string sLocalizerString = "LOCALIZER";
if (std::string(sTmpString).find(sLocalizerString) != std::string::npos) {
// std::cout << "Loading Localizer: "<<pcFName << '\n';
} else {
cerr << "ERROR: the image "<<pcFName <<" is not a LOCALIZER" << '\n';
free (sTmpString);
return -1;
}
/* read some useful tag values */
double
dIntensityWindow[2];
std::vector<std::string> tokens;
/* Intensity window center */
strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0028,0x1050), ds));
for (auto i = strtok(&sTmpString[0], "\\"); i != NULL; i = strtok(NULL, "\\"))
tokens.push_back(i);
if(tokens.size()>0){
dIntensityWindow[0] = std::atof(tokens.at(0).c_str() );
} else {
dIntensityWindow[0] = 0.;
}
tokens.clear();
/* Intensity window width */
strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0028,0x1051), ds));
for (auto i = strtok(&sTmpString[0], "\\"); i != NULL; i = strtok(NULL, "\\"))
tokens.push_back(i);
if(tokens.size()>0){
dIntensityWindow[1] = std::atof(tokens.at(0).c_str() );
} else {
dIntensityWindow[1] = 0.;
}
tokens.clear();
/* Image orientation */
eImageOrientationType
currImgOrient;
strcpy( sTmpString,gGetStringValueFromTag( gdcm::Tag(0x0018,0x5100), ds));
*std::remove(sTmpString, sTmpString + strlen(sTmpString), ' ') = 0;
if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::HFS])){
currImgOrient = eImageOrientationType::HFS;
} else if(!strcmp(sTmpString,ImageOrientationStrings[eImageOrientationType::FFS])){
currImgOrient = eImageOrientationType::FFS;
} else {
std::cerr<< "Image Orientation: unrecognised"<< sTmpString <<std::endl;
}
free (sTmpString);
/* Actual file read */
ImageReaderType2D::Pointer imageReader2D;
imageReader2D = ImageReaderType2D::New();
ImageIOType::Pointer dicomIO = ImageIOType::New();
imageReader2D->SetFileName(pcFName);
imageReader2D->SetImageIO(dicomIO);
try
{
imageReader2D->Update();
}
catch (const itk::ExceptionObject & ex)
{
std::cout << ex << std::endl;
return -1;
}
/* 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();
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];
break;
case eImageOrientationType::FFS:
toIECDirection.GetVnlMatrix()[iIdx][jIdx] = FFS2IEC[jIdx+iIdx*3];
break;
default:
break;
}
}
}
/* calculate image orientation with respect to fixed IEC */
InternalImageType::DirectionType LocalizerImagDirectionIEC =
toIECDirection * LocalizerImagDirectionDCM;
/* we calculate dot products between the Z components */
double dSumPA = 0;
dSumPA =
(LocalizerImagDirectionIEC[0][2] * PAElementsIEC[2]) +
(LocalizerImagDirectionIEC[1][2] * PAElementsIEC[5]) +
(LocalizerImagDirectionIEC[2][2] * PAElementsIEC[8]);
double dSumLAT = 0;
dSumLAT=
(LocalizerImagDirectionIEC[0][2] * LATElementsIEC[2]) +
(LocalizerImagDirectionIEC[1][2] * LATElementsIEC[5]) +
(LocalizerImagDirectionIEC[2][2] * LATElementsIEC[8]);
tProjOrientationType
currProjOrientation = NA;
/* 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"<<std::endl;
} else if (dSumLAT == 1 || dSumLAT == -1){
currProjOrientation = LAT;
// std::cout<<"It's a LAT image"<<std::endl;
} else {
return -1;
}
/* cast type to internal processing type */
CastFilterType3D::Pointer caster2DInput;
caster2DInput = CastFilterType3D::New();
caster2DInput->SetInput(imageReader2D->GetOutput());
caster2DInput->Update();
InternalImageType::Pointer MyLocalizerImage =
caster2DInput->GetOutput();
DuplicatorType::Pointer m_Duplicator;
TopogramImageMetaInformation::Pointer m_TImageMeta;
switch (currProjOrientation) {
case eProjectionOrientationType::PA:
if (m_PASourceDupli) {
m_PASourceDupli = NULL;
m_PASourceDupli = DuplicatorType::New();
}
m_Duplicator = m_PASourceDupli;
m_TImageMeta = m_TImage1MetaInfo;
break;
case eProjectionOrientationType::LAT:
if (m_LATSourceDupli) {
m_LATSourceDupli = NULL;
m_LATSourceDupli = DuplicatorType::New();
}
m_Duplicator = m_LATSourceDupli;
m_TImageMeta = m_TImage2MetaInfo;
break;
default:
return -1;
break;
}
if(m_TImageMeta != NULL){
//TODO
}
m_TImageMeta = TopogramImageMetaInformation::New();
m_TImageMeta->SetWLLevel(dIntensityWindow[0]);
m_TImageMeta->SetWLWindow(dIntensityWindow[1]);
m_TImageMeta->SetPatientOrientation(
currImgOrient);
m_TImageMeta->SetProjectionOrientation(
currProjOrientation);
m_Duplicator->SetInputImage(caster2DInput->GetOutput() );
m_Duplicator->Update();
return
(int) currProjOrientation;
}
double itkImageProcessor::GetLocalizerDisplayWindowLevel(int iImg)
{
TopogramImageMetaInformation::Pointer m_TImageMeta;
switch (iImg) {
case 0:
m_TImageMeta = m_TImage1MetaInfo;
break;
case 1:
m_TImageMeta = m_TImage2MetaInfo;
break;
default:
return 0;
break;
}
if(m_TImageMeta == NULL){
return 0;
} else {
return
m_TImageMeta->GetWLLevel();
}
return 0.;
}
double itkImageProcessor::GetLocalizerDisplayWindowWidth(int iImg)
{
TopogramImageMetaInformation::Pointer m_TImageMeta;
switch (iImg) {
case 0:
m_TImageMeta = m_TImage1MetaInfo;
break;
case 1:
m_TImageMeta = m_TImage2MetaInfo;
break;
default:
return 0;
break;
}
if(m_TImageMeta == NULL){
return 0;
} else {
return
m_TImageMeta->GetWLWindow();
}
return 0.;
}
Optimizer::ParametersType
itkImageProcessor::GetFinalR23Parameters(){
if (!m_TransformMetaInfo) {
itkExceptionMacro(<< "m_TransformMetaInfo not present");
}
Optimizer::ParametersType pPars(6);
pPars.fill(0.);
itk::Optimizer::ParametersType currentPosition =
m_R23->GetCurrentPosition();
switch (m_r23MetaInfo->GetDegreeOfFreedom()) {
case THREE_DOF:
pPars[3] = currentPosition[0]
- m_TransformMetaInfo->GetHiddenTranslations()[0];
pPars[4] = currentPosition[1]
- m_TransformMetaInfo->GetHiddenTranslations()[1];
pPars[5] = currentPosition[2]
- m_TransformMetaInfo->GetHiddenTranslations()[2];
break;
case SIX_DOF:
pPars[3] = currentPosition[0]
- m_TransformMetaInfo->GetHiddenTranslations()[0];
pPars[4] = currentPosition[1]
- m_TransformMetaInfo->GetHiddenTranslations()[1];
pPars[5] = currentPosition[2]
- m_TransformMetaInfo->GetHiddenTranslations()[2];
pPars[0] = currentPosition[3]
- m_TransformMetaInfo->GetHiddenRotations()[0];
pPars[1] = currentPosition[4]
- m_TransformMetaInfo->GetHiddenRotations()[1];
pPars[2] = currentPosition[5]
- m_TransformMetaInfo->GetHiddenRotations()[2];
break;
default:
pPars.fill(0.);
break;
}
return
pPars;
}
void itkImageProcessor::SetOptimizer(tOptimizerTypeEnum eOpti){
if(eOpti != tOptimizerTypeEnum::AMOEBA &&
eOpti != tOptimizerTypeEnum::POWELL){
itkExceptionMacro(<< "Unkown optimzer : " << eOpti);
return;
}
std::cout<<"Setting Optimizer: "<<eOpti<<std::endl;
m_r23MetaInfo->SetOptimizerType(eOpti);
}
void itkImageProcessor::SetMetric(tMetricTypeEnum eMetric){
if(eMetric != tMetricTypeEnum::NCC &&
eMetric != tMetricTypeEnum::MI){
itkExceptionMacro(<< "Unkown metric string : " << eMetric);
return;
}
std::cout<<"Setting Metric: "<<eMetric<<std::endl;
m_r23MetaInfo->SetMetricType(eMetric);
}
void itkImageProcessor::SetOptimizer(std::string optimizer)
{
if (optimizer.compare("Powell") == 0) {
m_r23MetaInfo->SetOptimizerType(tOptimizerTypeEnum::POWELL);
} else if (optimizer.compare("Amoeba") == 0) {
m_r23MetaInfo->SetOptimizerType(tOptimizerTypeEnum::AMOEBA);
} else if (optimizer.compare("Exhaustive") == 0) {
m_r23MetaInfo->SetOptimizerType(tOptimizerTypeEnum::EXHAUSTIVE);
} else {
itkExceptionMacro(<< "Unkown optimzer string : " << optimizer);
}
}
void itkImageProcessor::SetMetric(std::string metric)
{
if (metric.compare("NCC") == 0) {
m_r23MetaInfo->SetMetricType(tMetricTypeEnum::NCC);
} else if (metric.compare("MI") == 0) {
m_r23MetaInfo->SetMetricType(tMetricTypeEnum::MI);
} else {
itkExceptionMacro(<< "Unkown metric string : " << metric);
}
}
void itkImageProcessor::SetDegreeOfFreedom(tDegreeOfFreedomEnum dof){
m_DRTGeometryMetaInfo->SetDegreeOfFreedom(dof);
m_r23MetaInfo->SetDegreeOfFreedom(dof);
}
void itkImageProcessor::SetHandleRotationImportOffset(eHandlingRotImpTransform hrio){
m_DRTGeometryMetaInfo->SetHandleRotationImpOffset(hrio);
}
void itkImageProcessor::SetPowellOptimParameters(
double dStepT,
double dValTol,
double dStepL,
int iMaxLinI,
int iMaxIter){
m_PowellOptimizerMetaInfo->SetStepTolerance(dStepT);
m_PowellOptimizerMetaInfo->SetValueTolerance(dValTol);
m_PowellOptimizerMetaInfo->SetStepLength(dStepL);
m_PowellOptimizerMetaInfo->SetMaximumLineInteration(iMaxLinI);
m_PowellOptimizerMetaInfo->SetMaxIterations(iMaxIter);
}
void itkImageProcessor::SetAmoebaOptimParameters(
double dParConvTol,
double dFuntConvTol,
double dSimplex,
int iMaxIter){
m_AmoebaOptimizerMetaInfo->SetParametersConvergenceTolerance(dParConvTol);
m_AmoebaOptimizerMetaInfo->SetFunctionConvergenceTolerance(dFuntConvTol);
m_AmoebaOptimizerMetaInfo->SetSimplexDelta(dSimplex);
m_AmoebaOptimizerMetaInfo->SetMaxIterations(iMaxIter);
}
void itkImageProcessor::SetMIMetricParameters(double dMaxT,int iNhb){
m_MIMetricMetaInfo->SetMaxTranslation(dMaxT);
m_MIMetricMetaInfo->SetNumberOfHistogramBins(iNhb);
}
void itkImageProcessor::SetNCCMetricParameters(double dMaxT,bool bSm){
m_NCCMetricMetaInfo->SetMaxTranslation(dMaxT);
m_NCCMetricMetaInfo->SetSubtractMean(bSm);
}
void itkImageProcessor::InitializeProjector()
{
if(m_DRTImage1MetaInfo == NULL ||
m_DRTImage2MetaInfo == NULL ||
m_DRTGeometryMetaInfo == NULL ||
m_TransformMetaInfo == NULL ){
//TODO
return;
}
const double dtr = (atan(1.0) * 4.0) / 180.0;
ImageType3D::PointType ZeroPoint;
ZeroPoint.Fill(0.);
TransformType::Pointer CurrTransform;
CurrTransform= CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_InternalTransf1->GetpCalProjCenter(),
m_InternalTransf1->GetpRTIsocenter(),
m_InternalTransf1->GetIECtoLPSDirs());
transform1->SetComputeZYX(true);
transform1->SetIdentity();
transform1->SetTranslation(
CurrTransform->GetTranslation());
transform1->SetRotation(
CurrTransform->GetAngleX(),
CurrTransform->GetAngleY(),
CurrTransform->GetAngleZ()
);
transform1->SetCenter(
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero() );
// 2D Image 1
interpolator1->SetProjectionAngle(
dtr * m_DRTImage1MetaInfo->GetProjectionAngleLPS() );
interpolator1->SetFocalPointToIsocenterDistance(
m_DRTImage1MetaInfo->GetSCD());
interpolator1->SetThreshold(
m_DRTGeometryMetaInfo->GetIntensityThreshold() );
interpolator1->SetPanelOffset(
m_DRTImage1MetaInfo->GetPanelOffsetPGeo());
interpolator1->SetTransform(transform1);
interpolator1->Initialize();
CurrTransform= CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_InternalTransf2->GetpCalProjCenter(),
m_InternalTransf2->GetpRTIsocenter(),
m_InternalTransf2->GetIECtoLPSDirs());
transform2->SetComputeZYX(true);
transform2->SetIdentity();
transform2->SetTranslation(
CurrTransform->GetTranslation());
transform2->SetRotation(
CurrTransform->GetAngleX(),
CurrTransform->GetAngleY(),
CurrTransform->GetAngleZ()
);
transform2->SetCenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero() );
// 2D Image 2
interpolator2->SetProjectionAngle(
dtr * m_DRTImage2MetaInfo->GetProjectionAngleLPS() );
interpolator2->SetFocalPointToIsocenterDistance(
m_DRTImage2MetaInfo->GetSCD() );
interpolator2->SetThreshold(
m_DRTGeometryMetaInfo->GetIntensityThreshold() );
interpolator2->SetPanelOffset(
m_DRTImage2MetaInfo->GetPanelOffsetPGeo());
interpolator2->SetTransform(transform2);
interpolator2->Initialize();
resampleFilter1 = ResampleFilterType::New();
resampleFilter1->SetInput(
m_VolumeSourceDupli->GetOutput());
ImageType3D::RegionType lagerReg =
m_VolumeSourceDupli->GetOutput()->GetLargestPossibleRegion();
resampleFilter1->SetDefaultPixelValue(0);
resampleFilter1->SetNumberOfWorkUnits(iNWUnits);
// The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance
// have been set before registration. Here we only need to replace the initial
// transform with the final transform.
interpolator1->SetTransform(transform1);
interpolator1->Initialize();
resampleFilter1->SetInterpolator(interpolator1);
resampleFilter1->SetSize(m_DRTImage1MetaInfo->GetSize());
resampleFilter1->SetOutputSpacing(m_DRTImage1MetaInfo->GetSpacing());
resampleFilter1->SetOutputOrigin(m_DRTImage1MetaInfo->GetOrigin());
resampleFilter2 = ResampleFilterType::New();
resampleFilter2->SetInput(
m_VolumeSourceDupli->GetOutput());
resampleFilter2->SetDefaultPixelValue(0);
resampleFilter2->SetNumberOfWorkUnits(iNWUnits);
// The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance
// have been set before registration. Here we only need to replace the initial
// transform with the final transform.
interpolator2->SetTransform(transform2);
interpolator2->Initialize();
resampleFilter2->SetInterpolator(interpolator2);
resampleFilter2->SetSize(m_DRTImage2MetaInfo->GetSize());
resampleFilter2->SetOutputSpacing(m_DRTImage2MetaInfo->GetSpacing());
resampleFilter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOrigin());
filter1 =
ChangeInformationFilterType::New();
filter2 =
ChangeInformationFilterType::New();
/* First of all we need the center of the Projection images in its reference frame */
resampleFilter1->Update();
filter1->SetInput( resampleFilter1->GetOutput() );
filter1->SetOutputDirection(m_DRTImage1MetaInfo->GetImageDirectionsLPS() );//IECtoLPS_Directions);
filter1->ChangeDirectionOn();
filter1->SetOutputOrigin( m_DRTImage1MetaInfo->GetOriginLPS() );
filter1->ChangeOriginOn();
filter1->UpdateOutputInformation();
/* Again.. */
resampleFilter2->Update();
filter2->SetInput( resampleFilter2 ->GetOutput() );
filter2->SetOutputDirection(m_DRTImage2MetaInfo->GetImageDirectionsLPS() );
filter2->ChangeDirectionOn();
filter2->SetOutputOrigin(m_DRTImage2MetaInfo->GetOriginLPS() );
filter2->ChangeOriginOn();
filter2->UpdateOutputInformation();
filter1->Update();
filter2->Update();
imageDRT1In = filter1->GetOutput();
imageDRT2In = filter2->GetOutput();
}
itkImageProcessor::InternalImageType::DirectionType
itkImageProcessor::CalcDRTImageDirections(
double angle,
InternalImageType::DirectionType DRT2LPS){
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 *= DRT2LPS;
return
DRTDirs;
}
int itkImageProcessor::unloadRTPlanAndMeta(){
std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
m_RTMetaInfo = NULL;
ImageType3D::PointType pZero;
pZero.Fill(0.);
m_CTMetaInfo->SetImportOffset(pZero);
m_TransformMetaInfo->SetHiddenRotations(pZero);
return 0;
}
void itkImageProcessor::loadRTPlanInfo(
double dIsoX, double dIsoY, double dIsoZ,
double dLAT, double dVRT ,double dLNG){
if(m_RTMetaInfo != NULL){
std::cout << " NEVER HERE loadRTPlanInfo m_RTMetaInfo" << std::endl;
//TODO
}
if(m_CTMetaInfo != NULL){
std::cout << " ALWAYS HERE loadRTPlanInfo m_CTMetaInfo" << std::endl;
//TODO
}
if(m_DRTGeometryMetaInfo != NULL){
std::cout << " ALWAYS HERE loadRTPlanInfo m_DRTGeometryMetaInfo" << std::endl;
//TODO
}
m_RTMetaInfo = RTGeometryMetaInformation::New();
ImageType3D::PointType
Point;
Point[0] = dIsoX;
Point[1] = dIsoY;
Point[2] = dIsoZ;
m_RTMetaInfo->SetIsocenterLPS(Point);
Point[0] = dLAT;
Point[1] = dLNG;
Point[2] = dVRT;
m_RTMetaInfo->SetIsocenterIECS(Point);
ImageType3D::PointType
pZero (3);
pZero.Fill(0.);
m_CTMetaInfo->SetImportOffset(
CalcImportVolumeOffset(
m_RTMetaInfo->GetIsocenterIECS(),
m_CTMetaInfo->GetLPS2IECDirections(),
m_RTMetaInfo->GetIsocenterLPS(),
m_DRTGeometryMetaInfo->GetIECS2IECScannerT(),
(m_DRTGeometryMetaInfo->GetUseRotationsForHiddenTransform() ?
m_DRTGeometryMetaInfo->GetIECS2IECScannerR() : pZero)
)
);
m_TransformMetaInfo->SetHiddenRotations(
(m_DRTGeometryMetaInfo->GetUseRotationsForHiddenTransform() ?
m_DRTGeometryMetaInfo->GetIECS2IECScannerR() : pZero)
);
this->UpdateProjectionGeometryMeta();
this->InitializeProjector();
}
void itkImageProcessor::UpdateProjectionGeometryMeta(){
if(m_CTMetaInfo == NULL){
return;
//TODO.
}
if(m_DRTGeometryMetaInfo == NULL){
return;
//TODO.
}
if(m_DRTImage1MetaInfo == NULL){
return;
//TODO.
}
if(m_DRTImage2MetaInfo == NULL){
return;
//TODO.
}
// FIRST
ImageType3D::PointType NominalIsocenter_wZcorrectionLPS;
NominalIsocenter_wZcorrectionLPS =
m_CTMetaInfo->GetProjectionOriginLPS(
m_DRTGeometryMetaInfo->GetProjectionCenter()
);
ImageType3D::PointType NominalIsocenterZero_wZcorrectionLPS;
NominalIsocenterZero_wZcorrectionLPS =
m_CTMetaInfo->GetProjectionOriginLPSZero(
m_DRTGeometryMetaInfo->GetProjectionCenter()
);
ImageType3D::PointType IsocenterOffsetLPS;
IsocenterOffsetLPS = m_CTMetaInfo->ConvertIECPointToLPSPoint(
m_DRTGeometryMetaInfo->GetProjectionCenterOffset1());
ImageType3D::PointType CalibratedIsocenterLPS;
CalibratedIsocenterLPS[0] = NominalIsocenter_wZcorrectionLPS[0] +
IsocenterOffsetLPS[0];
CalibratedIsocenterLPS[1] = NominalIsocenter_wZcorrectionLPS[1] +
IsocenterOffsetLPS[1];
CalibratedIsocenterLPS[2] = NominalIsocenter_wZcorrectionLPS[2] +
IsocenterOffsetLPS[2];
ImageType3D::PointType CalibratedIsocenterZeroLPS;
CalibratedIsocenterZeroLPS[0] = NominalIsocenterZero_wZcorrectionLPS[0] +
IsocenterOffsetLPS[0] /*+ m_CTMetaInfo->GetSpacing()[0]/2*/;
CalibratedIsocenterZeroLPS[1] = NominalIsocenterZero_wZcorrectionLPS[1] +
IsocenterOffsetLPS[1] /*+ m_CTMetaInfo->GetSpacing()[1]/2*/;
CalibratedIsocenterZeroLPS[2] = NominalIsocenterZero_wZcorrectionLPS[2] +
IsocenterOffsetLPS[2] /*+ m_CTMetaInfo->GetSpacing()[2]/2*/;
m_DRTImage1MetaInfo->SetProjectionAngleLPS(
this->CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() +
m_DRTGeometryMetaInfo->GetProjectionAngle1OffsetIEC())
);
m_DRTImage1MetaInfo->SetSCD(
m_DRTGeometryMetaInfo->GetSCD1() +
m_DRTGeometryMetaInfo->GetSCD1Offset()
);
m_DRTImage1MetaInfo->SetProjectionOriginLPS(
CalibratedIsocenterLPS);
m_DRTImage1MetaInfo->SetProjectionOriginLPSZero(
CalibratedIsocenterZeroLPS);
m_DRTImage1MetaInfo->SetSizeWithBounds(
NULL,//vBounds.data(),
m_DRTGeometryMetaInfo->GetDRT1Size(),
m_DRTGeometryMetaInfo->GetDRT1Spacing()
);
m_DRTImage1MetaInfo->SetPanelOffsetPGeo(
this->CalcPanelOffsetPGeo(m_CTMetaInfo->GetPatientOrientation(),
-m_DRTGeometryMetaInfo->GetPanel1Offset())
);
// This HAS TO be calculated from the nominal isocenter.
m_DRTImage1MetaInfo->SetOriginLPS(
CalcDRTImageOrigin(
m_DRTImage1MetaInfo->GetOrigin(),
m_DRTImage1MetaInfo->GetCOV(),
NominalIsocenter_wZcorrectionLPS,
this->CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle1IEC() ),
Std_DRT2LPS
)
);
m_DRTImage1MetaInfo->SetImageDirectionsLPS(
CalcDRTImageDirections(
this->CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle1IEC()),
Std_DRT2LPS)
);
if(m_RTMetaInfo == NULL)
{
m_InternalTransf1->SetpCalProjCenter(
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero());
m_InternalTransf1->SetpRTIsocenter(
m_DRTImage1MetaInfo->GetProjectionOriginLPS());
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions =
m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
m_InternalTransf1->SetIECtoLPSDirs(IECtoLPS_Directions);
} else {
m_InternalTransf1->SetpCalProjCenter(
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero());
m_InternalTransf1->SetpRTIsocenter(
m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS());
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions =
m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
m_InternalTransf1->SetIECtoLPSDirs(IECtoLPS_Directions);
}
// SECOND
IsocenterOffsetLPS = m_CTMetaInfo->ConvertIECPointToLPSPoint(
m_DRTGeometryMetaInfo->GetProjectionCenterOffset2());
CalibratedIsocenterLPS[0] = NominalIsocenter_wZcorrectionLPS[0] +
IsocenterOffsetLPS[0];
CalibratedIsocenterLPS[1] = NominalIsocenter_wZcorrectionLPS[1] +
IsocenterOffsetLPS[1];
CalibratedIsocenterLPS[2] = NominalIsocenter_wZcorrectionLPS[2] +
IsocenterOffsetLPS[2];
CalibratedIsocenterZeroLPS[0] = NominalIsocenterZero_wZcorrectionLPS[0] +
IsocenterOffsetLPS[0];
CalibratedIsocenterZeroLPS[1] = NominalIsocenterZero_wZcorrectionLPS[1] +
IsocenterOffsetLPS[1];
CalibratedIsocenterZeroLPS[2] = NominalIsocenterZero_wZcorrectionLPS[2] +
IsocenterOffsetLPS[2];
m_DRTImage2MetaInfo->SetProjectionAngleLPS(
this->CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle2IEC() +
m_DRTGeometryMetaInfo->GetProjectionAngle2OffsetIEC())
);
m_DRTImage2MetaInfo->SetSCD(
m_DRTGeometryMetaInfo->GetSCD2()+
m_DRTGeometryMetaInfo->GetSCD2Offset()
);
m_DRTImage2MetaInfo->SetProjectionOriginLPS(
CalibratedIsocenterLPS);
m_DRTImage2MetaInfo->SetProjectionOriginLPSZero(
CalibratedIsocenterZeroLPS);
m_DRTImage2MetaInfo->SetSizeWithBounds(
NULL,//vBounds.data(),
m_DRTGeometryMetaInfo->GetDRT2Size(),
m_DRTGeometryMetaInfo->GetDRT2Spacing()
);
m_DRTImage2MetaInfo->SetPanelOffsetPGeo(
this->CalcPanelOffsetPGeo(m_CTMetaInfo->GetPatientOrientation(),
-m_DRTGeometryMetaInfo->GetPanel2Offset()));
m_DRTImage2MetaInfo->SetOriginLPS(
CalcDRTImageOrigin(
m_DRTImage2MetaInfo->GetOrigin(),
m_DRTImage2MetaInfo->GetCOV(),
NominalIsocenter_wZcorrectionLPS,
this->CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) ,
Std_DRT2LPS
));
m_DRTImage2MetaInfo->SetImageDirectionsLPS(
CalcDRTImageDirections(
this->CalcProjectionAngleLPS(
m_CTMetaInfo->GetPatientOrientation(),
m_DRTGeometryMetaInfo->GetProjectionAngle2IEC()) ,
Std_DRT2LPS)
);
if(m_RTMetaInfo == NULL)
{
m_InternalTransf2->SetpCalProjCenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero());
m_InternalTransf2->SetpRTIsocenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPS());
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions =
m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
m_InternalTransf2->SetIECtoLPSDirs(IECtoLPS_Directions);
} else {
m_InternalTransf2->SetpCalProjCenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero());
m_InternalTransf2->SetpRTIsocenter(
m_RTMetaInfo->GetIsocenterLPS() - m_CTMetaInfo->GetOriginLPS());
InternalImageType::DirectionType IECtoLPS_Directions;
IECtoLPS_Directions =
m_CTMetaInfo->GetLPS2IECDirections().GetTranspose();
m_InternalTransf2->SetIECtoLPSDirs(IECtoLPS_Directions);
}
}
/*
Nobody should go through all this...
To be very careful editing ...
*/
itkImageProcessor::ImageType3D::PointType
itkImageProcessor::CalcDRTImageOrigin(
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);
ImageType3D::PointType NewOrigin =
m_ProjectionCenterLPS + DRT2LPS * (m_DRTOrigin - m_DRTCOV);
return
NewOrigin;
}
void itkImageProcessor::GetProjectionImages(){
if(m_DRTImage1MetaInfo == NULL ||
m_DRTImage2MetaInfo == NULL ||
m_DRTGeometryMetaInfo == NULL ||
m_TransformMetaInfo == NULL ){
return;
//TODO
return;
}
ImageType3D::PointType ZeroPoint;
ZeroPoint.Fill(0.);
transform1->SetComputeZYX(true);
transform1->SetIdentity();
TransformType::Pointer CurrTransform;
CurrTransform= CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_InternalTransf1->GetpCalProjCenter(),
m_InternalTransf1->GetpRTIsocenter(),
m_InternalTransf1->GetIECtoLPSDirs());
transform1->SetComputeZYX(true);
transform1->SetIdentity();
transform1->SetTranslation(
CurrTransform->GetTranslation());
transform1->SetRotation(
CurrTransform->GetAngleX(),
CurrTransform->GetAngleY(),
CurrTransform->GetAngleZ()
);
transform1 ->SetCenter(
m_DRTImage1MetaInfo->GetProjectionOriginLPSZero());
// The parameters of interpolator1, such as ProjectionAngle and FocalPointToIsocenterDistance
// have been set before registration. Here we only need to replace the initial
// transform with the final transform.
interpolator1->SetTransform(transform1);
interpolator1->Initialize();
transform2->SetComputeZYX(true);
transform2->SetIdentity();
CurrTransform = CalculateInternalTransformV3(
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR(),
m_InternalTransf2->GetpCalProjCenter(),
m_InternalTransf2->GetpRTIsocenter(),
m_InternalTransf2->GetIECtoLPSDirs());
transform2->SetComputeZYX(true);
transform2->SetIdentity();
transform2->SetTranslation(
CurrTransform->GetTranslation());
transform2->SetRotation(
CurrTransform->GetAngleX(),
CurrTransform->GetAngleY(),
CurrTransform->GetAngleZ()
);
transform2 ->SetCenter(
m_DRTImage2MetaInfo->GetProjectionOriginLPSZero());
// // The parameters of interpolator2, such as ProjectionAngle and FocalPointToIsocenterDistance
// // have been set before registration. Here we only need to replace the initial
// // transform with the final transform.
interpolator2->SetTransform(transform2); //finalTransform);
interpolator2->Initialize();
filter1->Update();
filter2->Update();
imageDRT1In = filter1->GetOutput();
imageDRT2In = filter2->GetOutput();
}
itkImageProcessor::ImageType3D::PointType
itkImageProcessor::CalcImportVolumeOffset(
ImageType3D::PointType rtCouchOffset,
InternalImageType::DirectionType VolumeLPStoIEC,
ImageType3D::PointType rtIsocenterLPS,
ImageType3D::PointType IEC2DCMMapT,
ImageType3D::PointType IEC2DCMMapR){
TransformType::Pointer InputTransform = TransformType::New();
InputTransform->SetComputeZYX(true);
InputTransform->SetIdentity();
TransformType::OutputVectorType translation;
translation[0] = IEC2DCMMapT[0];
translation[1] = IEC2DCMMapT[1];
translation[2] = IEC2DCMMapT[2];
InputTransform->SetTranslation(translation);
const double dtr = (atan(1.0) * 4.0) / 180.0;
InputTransform->SetRotation(
dtr * IEC2DCMMapR[0],
dtr * IEC2DCMMapR[1],
dtr * IEC2DCMMapR[2]);
ImageType3D::PointType m_TransformOrigin;
m_TransformOrigin.Fill(0.);
InputTransform->SetCenter(
m_TransformOrigin );
ImageType3D::PointType IsoSupport_IECPos =
InputTransform->TransformPoint(rtCouchOffset );
InternalImageType::DirectionType VolumeIECtoLPS;
VolumeIECtoLPS = VolumeLPStoIEC.GetTranspose();
ImageType3D::PointType IsoSupport_LPSPos =
VolumeIECtoLPS * IsoSupport_IECPos;
ImageType3D::PointType Offset =
rtIsocenterLPS - IsoSupport_LPSPos;
std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ <<
" " << Offset << std::endl;
return
Offset;
}
void itkImageProcessor::Write2DImages(){
// using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
// using WriterType = itk::ImageFileWriter<OutputImageType>;
// 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());
// WriterType::Pointer writer1 = WriterType::New();
// 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;
// }
// }
// 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();
// 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;
// }
// }
}
vtkImageData* itkImageProcessor::GetLocalizer1VTK()
{
using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
rescaler1->SetOutputMinimum(0);
rescaler1->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
rescaler1->SetInput( m_PASourceDupli->GetOutput() );
rescaler1->Update();
toVTKLocalizer1->SetInput(rescaler1->GetOutput());
toVTKLocalizer1->Update();
return
toVTKLocalizer1->GetOutput();
}
vtkImageData* itkImageProcessor::GetLocalizer2VTK()
{
using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
rescaler2->SetOutputMinimum(0);
rescaler2->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
rescaler2->SetInput( m_LATSourceDupli->GetOutput() );
rescaler2->Update();
toVTKLocalizer2->SetInput(rescaler2->GetOutput());
toVTKLocalizer2->Update();
return
toVTKLocalizer2->GetOutput();
}
vtkImageData* itkImageProcessor::GetProjection1VTK()
{
if(m_DRTImage1MetaInfo == NULL ||
m_DRTGeometryMetaInfo == NULL ||
m_TransformMetaInfo == NULL ){
return NULL;
//TODO
}
using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
rescaler1->SetOutputMinimum(0);
rescaler1->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
rescaler1->SetInput( imageDRT1In );
rescaler1->Update();
toVTK2D1->SetInput(rescaler1->GetOutput());
toVTK2D1->Update();
return
toVTK2D1->GetOutput();
}
vtkImageData* itkImageProcessor::GetProjection1VTKToWrite()
{
if(m_DRTImage1MetaInfo == NULL ||
m_DRTGeometryMetaInfo == NULL ||
m_TransformMetaInfo == NULL ){
return NULL;
//TODO
}
using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<InternalImageType>;
auto imageCalculatorFilter = ImageCalculatorFilterType::New();
imageCalculatorFilter->SetImage(imageDRT1In);
imageCalculatorFilter->Compute();
using IntWindowType = itk::IntensityWindowingImageFilter<InternalImageType, OutputImageType>;
auto intWindowFilter = IntWindowType::New();
intWindowFilter->SetWindowMinimum(imageCalculatorFilter->GetMinimum());
intWindowFilter->SetWindowMaximum(imageCalculatorFilter->GetMaximum());
intWindowFilter->SetOutputMinimum(0);
if(imageCalculatorFilter->GetMaximum() > SHRT_MAX){
intWindowFilter->SetOutputMaximum(SHRT_MAX);
}else
intWindowFilter->SetOutputMaximum(imageCalculatorFilter->GetMaximum());
intWindowFilter->SetInput(imageDRT1In);
intWindowFilter->Update();
toVTK2D1->SetInput(intWindowFilter->GetOutput());
toVTK2D1->Update();
return
toVTK2D1->GetOutput();
}
vtkMatrix4x4 * itkImageProcessor::GetProjection1VTKTransform()
{
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]);
}
}
m_Projection1VTKTransform->SetElement(0,3,
interpolator1->GetInverseTransform()->GetTranslation()[0]);
m_Projection1VTKTransform->SetElement(1,3,
interpolator1->GetInverseTransform()->GetTranslation()[1]);
m_Projection1VTKTransform->SetElement(2,3,
interpolator1->GetInverseTransform()->GetTranslation()[2]);
return
m_Projection1VTKTransform;
}
vtkMatrix4x4 * itkImageProcessor::GetProjection2VTKTransform()
{
m_Projection2VTKTransform->Identity();
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]);
}
}
m_Projection2VTKTransform->SetElement(0,3,
interpolator2->GetInverseTransform()->GetTranslation()[0]);
m_Projection2VTKTransform->SetElement(1,3,
interpolator2->GetInverseTransform()->GetTranslation()[1]);
m_Projection2VTKTransform->SetElement(2,3,
interpolator2->GetInverseTransform()->GetTranslation()[2]);
return
m_Projection2VTKTransform;
}
vtkImageData* itkImageProcessor::GetProjection2VTK()
{
if(m_DRTImage2MetaInfo == NULL ||
m_DRTGeometryMetaInfo == NULL ||
m_TransformMetaInfo == NULL ){
return NULL;
//TODO
}
using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
rescaler2->SetOutputMinimum(0);
rescaler2->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
rescaler2->SetInput( imageDRT2In );
rescaler2->Update();
toVTK2D2->SetInput(rescaler2->GetOutput());
toVTK2D2->Update();
return
toVTK2D2->GetOutput();
}
vtkImageData* itkImageProcessor::GetProjection2VTKToWrite()
{
if(m_DRTImage2MetaInfo == NULL ||
m_DRTGeometryMetaInfo == NULL ||
m_TransformMetaInfo == NULL ){
return NULL;
//TODO
}
using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<InternalImageType>;
auto imageCalculatorFilter = ImageCalculatorFilterType::New();
imageCalculatorFilter->SetImage(imageDRT2In);
imageCalculatorFilter->Compute();
using IntWindowType = itk::IntensityWindowingImageFilter<InternalImageType, OutputImageType>;
auto intWindowFilter = IntWindowType::New();
intWindowFilter->SetWindowMinimum(imageCalculatorFilter->GetMinimum());
intWindowFilter->SetWindowMaximum(imageCalculatorFilter->GetMaximum());
intWindowFilter->SetOutputMinimum(0);
if(imageCalculatorFilter->GetMaximum() > SHRT_MAX)
intWindowFilter->SetOutputMaximum(SHRT_MAX);
else
intWindowFilter->SetOutputMaximum(imageCalculatorFilter->GetMaximum());
intWindowFilter->SetInput(imageDRT2In);
intWindowFilter->Update();
toVTK2D2->SetInput(intWindowFilter->GetOutput());
toVTK2D2->Update();
return
toVTK2D2->GetOutput();
}
void itkImageProcessor::WriteProjectionImages()
{
// Rescale the intensity of the projection images to 0-255 for output.
using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
rescaler1->SetOutputMinimum(0);
rescaler1->SetOutputMaximum(255);
rescaler1->SetInput( imageDRT1In );
RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
rescaler2->SetOutputMinimum(0);
rescaler2->SetOutputMaximum(255);
rescaler2->SetInput( imageDRT2In );
using WriterType = itk::ImageFileWriter<OutputImageType>;
WriterType::Pointer writer1 = WriterType::New();
WriterType::Pointer writer2 = WriterType::New();
writer1->SetFileName("/mnt/gdrive/Scratch/gfattori/Projection1.tif");
writer1->SetInput(rescaler1->GetOutput());
try
{
writer1->Update();
}
catch (itk::ExceptionObject & err)
{
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
writer2->SetFileName("/mnt/gdrive/Scratch/gfattori/Projection2.tif");
writer2->SetInput(rescaler2->GetOutput());
try
{
writer2->Update();
}
catch (itk::ExceptionObject & err)
{
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
}
void itkImageProcessor::SetInitialTranslations(double dX, double dY, double dZ)
{
if(m_TransformMetaInfo == NULL){
//todo
}
ImageType3D::PointType Translations;
Translations[0]= dX;
Translations[1]= dY;
Translations[2]= dZ;
m_TransformMetaInfo->SetUserTranslations(Translations);
}
void itkImageProcessor::SetInitialRotations(double dX, double dY, double dZ)
{
ImageType3D::PointType Rotations;
Rotations[0]= dX;
Rotations[1]= dY;
Rotations[2]= dZ;
m_TransformMetaInfo->SetUserRotations(Rotations);
}
Optimizer::ParametersType
itkImageProcessor::GetUserIsocentricTransform(){
Optimizer::ParametersType
m_Pars(6);
m_Pars.fill(0.);
if(!m_TransformMetaInfo){
return m_Pars;
}
m_Pars[0] = m_TransformMetaInfo->GetUserRotations()[0];
m_Pars[1] = m_TransformMetaInfo->GetUserRotations()[1];
m_Pars[2] = m_TransformMetaInfo->GetUserRotations()[2];
m_Pars[3] = m_TransformMetaInfo->GetUserTranslations()[0];
m_Pars[4] = m_TransformMetaInfo->GetUserTranslations()[1];
m_Pars[5] = m_TransformMetaInfo->GetUserTranslations()[2];
return
m_Pars;
}
TransformType::Pointer
itkImageProcessor::GetCompleteIsocentricTransform(){
if(m_TransformMetaInfo == nullptr ||
m_CTMetaInfo == nullptr){
return nullptr;
}
ImageType3D::PointType mISORTIEC =
m_CTMetaInfo->GetLPS2IECDirections() * m_RTMetaInfo->GetIsocenterLPS();
mISORTIEC[0] *= -1;
mISORTIEC[1] *= -1;
mISORTIEC[2] *= -1;
TransformType::Pointer mTransf=
MapTransformToNewOrigin(
mISORTIEC,
m_TransformMetaInfo->GetT(),
m_TransformMetaInfo->GetR()
);
itk::TransformMetaInformation::PointType
pTM = m_CTMetaInfo->GetLPS2IECDirections() * m_CTMetaInfo->GetImportOffset();
TransformType::OutputVectorType pTransl;
pTransl[0] = mTransf->GetTranslation()[0] - pTM[0];
pTransl[1] = mTransf->GetTranslation()[1] - pTM[1];
pTransl[2] = mTransf->GetTranslation()[2] - pTM[2];
mTransf->SetTranslation(pTransl);
itk::Point<double,3>pZero;
pZero.Fill(0.);
mTransf->SetCenter(pZero);
return mTransf;
}
void itkImageProcessor::SetRegionFixed1(
int iIdx0, int iIdx1, int iSz0, int iSz1){
auto region1temp = m_PASourceDupli->GetOutput()->GetBufferedRegion();
auto index1 = region1temp.GetIndex();
auto size1 = region1temp.GetSize();
index1[0] = iIdx0;
index1[1] = iIdx1;
size1[0] = iSz0;
size1[1] = iSz1;
size1[2] = 1;
roiAutoReg1.SetIndex(index1);
roiAutoReg1.SetSize(size1);
this->m_R23->SetroiAutoReg1(roiAutoReg1);
}
void itkImageProcessor::SetRegionFixed2(
int iIdx0, int iIdx1, int iSz0, int iSz1){
auto region2temp = m_LATSourceDupli->GetOutput()->GetBufferedRegion();
auto index2 = region2temp.GetIndex();
auto size2 = region2temp.GetSize();
index2[0] = iIdx0;
index2[1] = iIdx1;
size2[0] = iSz0;
size2[1] = iSz1;
size2[2] = 1;
roiAutoReg2.SetIndex(index2);
roiAutoReg2.SetSize(size2);
this->m_R23->SetroiAutoReg2(roiAutoReg2);
}
void itkImageProcessor::ResetROIRegions(){
this->m_R23->setROISizeToZero();
}
void itkImageProcessor::InizializeRegistration(){
m_R23->SetPA(m_PASourceDupli->GetOutput());
m_R23->SetLAT(m_LATSourceDupli->GetOutput());
m_R23->SetVolume(m_VolumeSourceDupli->GetOutput());
m_R23->Setr23Meta(m_r23MetaInfo);
m_R23->SetPowellMeta(m_PowellOptimizerMetaInfo);
m_R23->SetAmoebaMeta(m_AmoebaOptimizerMetaInfo);
m_R23->SetMIMeta(m_MIMetricMetaInfo);
m_R23->SetNCCMeta(m_NCCMetricMetaInfo);
m_R23->SetInternalTransf1(m_InternalTransf1);
m_R23->SetInternalTransf2(m_InternalTransf2);
m_R23->Setfilter1(filter1);
m_R23->Setfilter2(filter2);
m_R23->Setinterpolator1(interpolator1);
m_R23->Setinterpolator2(interpolator2);
m_R23->SetTransformMetaInfo(m_TransformMetaInfo);
m_R23->InitializeRegistration();
}
} // end namespace itk
Reference in New Issue View Git Blame Copy Permalink