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

Merge branch 'ScoutRTUIDevel' into 'master'

Browse Source 
Pre-release for devs

Closes #23, #33, #36, #50, #47, and #46

See merge request cpt_bioeng/drt!30
This commit is contained in:
fattori_g
2023-05-23 17:07:47 +02:00
parent d05205ca88 205f2f9db1
commit 14f4aca6d1
24 changed files with 4714 additions and 946 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
reg23Topograms/itkDTRrecon/CMakeLists.txt
View File

@ -2,6 +2,8 @@ cmake_minimum_required(VERSION 2.8.12)
SET(LIB_NAME itkDTRrecon)
add_subdirectory(autoreg)
find_package(ITK REQUIRED)
include(${ITK_USE_FILE})
INCLUDE_DIRECTORIES(${ITK_INCLUDE_DIRS})
@ -14,17 +16,21 @@ INCLUDE_DIRECTORIES(${GDCM_INCLUDE_DIRS})
SET(SRCS
itkImageProcessor.cpp
itkImageProcessorHelpers.cpp
vtkContourTopogramProjectionFilter.cxx
DRTMetaInformation.cpp
itkReg23.cpp
)
SET(HDR
itkImageProcessor.h
itkImageProcessorHelpers.h
itkQtIterationUpdate.h
itkgSiddonJacobsRayCastInterpolateImageFunction.h
itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
vtkContourTopogramProjectionFilter.h
DRTMetaInformation.h
itkReg23.h
)
ADD_LIBRARY(${LIB_NAME} ${SRCS} ${HDR})
@ -34,6 +40,7 @@ SET(LINK_LIBS
${VTK_LIBRARIES}
${ITK_LIBRARIES}
gdcmMSFF
autoreg
)
TARGET_LINK_LIBRARIES(

187
reg23Topograms/itkDTRrecon/DRTMetaInformation.cpp
View File

@ -137,69 +137,67 @@ void DRTImageMetaInformation::SetSizeWithBounds(
SpacingType Spacing)
{
double dXmin = dBounds[0];
double dYmin = dBounds[1];
double dZmin = dBounds[2];
double dXmax = dBounds[3];
double dYmax = dBounds[4];
double dZmax = dBounds[5];
// double dXmin = dBounds[0];
// double dYmin = dBounds[1];
// double dZmin = dBounds[2];
// double dXmax = dBounds[3];
// double dYmax = dBounds[4];
// double dZmax = dBounds[5];
double
dNPixMin1 = 0.,
dNPixMin2 = 0.,
dNPixMax1 = 0.,
dNPixMax2 = 0.;
// double
// dNPixMin1 = 0.,
// dNPixMin2 = 0.,
// dNPixMax1 = 0.,
// dNPixMax2 = 0.;
dNPixMin1 = abs(dXmin / Spacing[0]);
dNPixMin2 = abs(dYmin / Spacing[0]);
dNPixMax1 = abs(dXmax / Spacing[0]);
dNPixMax2 = abs(dYmax / Spacing[0]);
// dNPixMin1 = abs(dXmin / Spacing[0]);
// dNPixMin2 = abs(dYmin / Spacing[0]);
// dNPixMax1 = abs(dXmax / Spacing[0]);
// dNPixMax2 = abs(dYmax / Spacing[0]);
double dPixMaxOverlap =dNPixMin1;
// double dPixMaxOverlap =dNPixMin1;
if(dPixMaxOverlap > dNPixMin2){
dPixMaxOverlap = dNPixMin2;
}
// if(dPixMaxOverlap > dNPixMin2){
// dPixMaxOverlap = dNPixMin2;
// }
if(dPixMaxOverlap > dNPixMax1){
dPixMaxOverlap = dNPixMax1;
}
// if(dPixMaxOverlap > dNPixMax1){
// dPixMaxOverlap = dNPixMax1;
// }
if(dPixMaxOverlap > dNPixMax2){
dPixMaxOverlap = dNPixMax2;
}
// if(dPixMaxOverlap > dNPixMax2){
// dPixMaxOverlap = dNPixMax2;
// }
// if(floor(dPixMaxOverlap*2) > MaxSize[0]) {
// m_Size[0] = MaxSize[0];
// } else {
// m_Size[0] = floor(dPixMaxOverlap*2);
// }
// dNPixMin1 = abs(dZmin / Spacing[1]);
// dNPixMax1 = abs(dZmax / Spacing[1]);
// dPixMaxOverlap =dNPixMin1;
// if(dPixMaxOverlap > dNPixMax1){
// dPixMaxOverlap = dNPixMax1;
// }
// if(floor(dPixMaxOverlap*2) > MaxSize[1]) {
// m_Size[1] = MaxSize[1];
// } else {
// m_Size[1] = floor(dPixMaxOverlap*2);
// }
if(floor(dPixMaxOverlap*2) > MaxSize[0]) {
m_Size[0] = MaxSize[0];
} else {
m_Size[0] = floor(dPixMaxOverlap*2);
}
dNPixMin1 = abs(dZmin / Spacing[1]);
dNPixMax1 = abs(dZmax / Spacing[1]);
dPixMaxOverlap =dNPixMin1;
if(dPixMaxOverlap > dNPixMax1){
dPixMaxOverlap = dNPixMax1;
}
if(floor(dPixMaxOverlap*2) > MaxSize[1]) {
m_Size[1] = MaxSize[1];
} else {
m_Size[1] = floor(dPixMaxOverlap*2);
}
m_Size[2] = 1;
m_Spacing = Spacing;
return;
}
@ -452,6 +450,10 @@ DRTProjectionGeometryImageMetaInformation(){
this->m_SCD2Offset = 0.;
this->m_Panel1Offset = 0.;
this->m_Panel2Offset = 0.;
this->m_IntensityThreshold=0.;
this->m_DRT1Size.Fill(0.);
@ -472,6 +474,8 @@ DRTProjectionGeometryImageMetaInformation(){
this->m_ProjectionCenter.Fill(0.);
this->m_UseRotationsForHiddenTransform = true;
}
void
@ -515,23 +519,100 @@ RTGeometryMetaInformation
}
InternalTransformMetaInformation::InternalTransformMetaInformation(){
m_pCalProjCenter.Fill(0.);
m_pRTIsocenter.Fill(0.);
m_IECtoLPSDirs.SetIdentity();
}
void
InternalTransformMetaInformation
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
InternalTransformMetaInformation
::~InternalTransformMetaInformation ()
{
}
R23MetaInformation::
R23MetaInformation (){
m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN;
m_OptimizerType = tOptimizerTypeEnum::POWELL;
m_MetricType = tMetricTypeEnum::NCC;
m_MaxIterations = 6;
}
void
R23MetaInformation
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
R23MetaInformation
::~R23MetaInformation ()
{
}
TransformMetaInformation ::
TransformMetaInformation (){
m_Translations.Fill(0.);
m_HiddenTranslations.Fill(0.);
m_Rotations.Fill(0.);
m_HiddenRotations.Fill(0.);
m_UserRotations.Fill(0.);
m_UserTranslations.Fill(0.);
m_ReferenceTransform.Fill(0.);
//m_ReferenceTransform.Fill(0.);
m_CurrentTransform.Fill(0.);
//m_CurrentTransform.Fill(0.);
// m_DegreeOfFreedom = tDegreeOfFreedomEnum::UNKNOWN;
}
TransformMetaInformation::PointType
TransformMetaInformation::GetT(){
PointType
pT;
pT[0] = m_UserTranslations[0] + m_HiddenTranslations[0];
pT[1] = m_UserTranslations[1] + m_HiddenTranslations[1];
pT[2] = m_UserTranslations[2] + m_HiddenTranslations[2];
return
pT;
}
TransformMetaInformation::PointType
TransformMetaInformation::TransformMetaInformation::GetR(){
PointType
pR;
pR[0] = m_UserRotations[0] + m_HiddenRotations[0];
pR[1] = m_UserRotations[1] + m_HiddenRotations[1];
pR[2] = m_UserRotations[2] + m_HiddenRotations[2];
return
pR;
}

239
reg23Topograms/itkDTRrecon/DRTMetaInformation.h
View File

@ -11,15 +11,83 @@
namespace itk
{
typedef enum eProjectionOrientationType
{NA = 0, LAT=2, PA=1}
tProjOrientationType;
typedef enum eProjectionOrientationType{
NA = 0,
LAT=2,
PA=1
} tProjOrientationType;
/** Enum type for Orientation */
typedef enum eImageOrientationType
{NotDefined = 0, HFS = 1, FFS = 2}
tPatOrientation;
typedef enum eImageOrientationType{
NotDefined = 0,
HFS = 1,
FFS = 2
} tPatOrientation;
typedef enum eDegreeOfFreedomType {
UNKNOWN = 0,
THREE_DOF,
SIX_DOF,
X_ONLY,
Y_ONLY,
Z_ONLY,
ROTATION_ONLY,
OTHER
} tDegreeOfFreedomEnum;
typedef enum eOptimizerType{
POWELL,
AMOEBA,
EXHAUSTIVE
} tOptimizerTypeEnum;
typedef enum eMetricType{
NCC,
MI
} tMetricTypeEnum;
inline int GetNumberOfDOF(eDegreeOfFreedomType dof)
{
switch (dof) {
case SIX_DOF:
return 6;
break;
case THREE_DOF:
case ROTATION_ONLY:
return 3;
break;
case X_ONLY:
case Y_ONLY:
case Z_ONLY:
return 1;
break;
default:
return 0;
}
}
class DegreeOfFreedom : public itk::Object {
public:
typedef DegreeOfFreedom Self;
typedef itk::Object Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro(Self);
protected:
bool
m_TranslationX,
m_TranslationY,
m_TranslationZ,
m_RotationX,
m_RotationY,
m_RotationZ;
tDegreeOfFreedomEnum m_dof;
};
class TopogramImageMetaInformation :
public itk::Object{
@ -317,6 +385,12 @@ public:
itkSetMacro(SCD2Offset, double);
itkGetMacro(SCD2Offset, double);
itkSetMacro(Panel1Offset, double);
itkGetMacro(Panel1Offset, double);
itkSetMacro(Panel2Offset, double);
itkGetMacro(Panel2Offset, double);
itkSetMacro(IntensityThreshold, double);
itkGetMacro(IntensityThreshold, double);
@ -347,6 +421,10 @@ public:
itkSetMacro(ProjectionCenter, PointType);
itkGetMacro(ProjectionCenter, PointType);
itkSetMacro(UseRotationsForHiddenTransform, bool);
itkGetMacro(UseRotationsForHiddenTransform, bool);
protected:
double
@ -358,6 +436,8 @@ protected:
m_SCD2,
m_SCD1Offset,
m_SCD2Offset,
m_Panel1Offset,
m_Panel2Offset,
m_IntensityThreshold;
@ -384,6 +464,9 @@ protected:
m_ProjectionCenterOffset2;
bool
m_UseRotationsForHiddenTransform;
/** Default Constructor **/
DRTProjectionGeometryImageMetaInformation ();
/** Default Destructor **/
@ -451,6 +534,126 @@ private:
};
class InternalTransformMetaInformation :
public itk::Object{
public:
/** standard typedefs **/
typedef InternalTransformMetaInformation Self;
typedef itk::Object Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::Point<double, 3> PointType;
typedef itk::Matrix<double,3,3> TransformMatrixType;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(InternalTransformMetaInformation, itk::Object);
/** object information streaming **/
void PrintSelf(std::ostream& os, itk::Indent indent) const;
itkSetMacro(pCalProjCenter,PointType);
itkSetMacro(pRTIsocenter,PointType);
itkSetMacro(IECtoLPSDirs,TransformMatrixType);
itkGetMacro(pCalProjCenter,PointType);
itkGetMacro(pRTIsocenter,PointType);
itkGetMacro(IECtoLPSDirs,TransformMatrixType);
protected:
PointType
m_pCalProjCenter,
m_pRTIsocenter;
TransformMatrixType
m_IECtoLPSDirs;
/** Default Constructor **/
InternalTransformMetaInformation ();
/** Default Destructor **/
virtual ~InternalTransformMetaInformation ();
private:
/** purposely not implemented **/
InternalTransformMetaInformation (const Self&);
/** purposely not implemented **/
void operator=(const Self&);
};
class R23MetaInformation :
public itk::Object{
public:
/** standard typedefs **/
typedef R23MetaInformation Self;
typedef itk::Object Superclass;
typedef itk::SmartPointer<Self> Pointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(R23MetaInformation, itk::Object);
/** object information streaming **/
void PrintSelf(std::ostream& os, itk::Indent indent) const override;
itkSetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum);
itkGetEnumMacro(DegreeOfFreedom, tDegreeOfFreedomEnum);
itkSetEnumMacro(OptimizerType, tOptimizerTypeEnum);
itkGetEnumMacro(OptimizerType, tOptimizerTypeEnum);
itkSetEnumMacro(MetricType, tMetricTypeEnum);
itkGetEnumMacro(MetricType, tMetricTypeEnum);
itkSetEnumMacro(MaxIterations, int);
itkGetEnumMacro(MaxIterations, int);
protected:
tDegreeOfFreedomEnum
m_DegreeOfFreedom;
tOptimizerTypeEnum
m_OptimizerType;
tMetricTypeEnum
m_MetricType;
int
m_MaxIterations;
/** Default Constructor **/
R23MetaInformation ();
/** Default Destructor **/
virtual ~R23MetaInformation ();
private:
/** purposely not implemented **/
R23MetaInformation (const Self&);
/** purposely not implemented **/
void operator=(const Self&);
};
class TransformMetaInformation :
public itk::Object{
@ -472,11 +675,11 @@ public:
void PrintSelf(std::ostream& os, itk::Indent indent) const;
itkSetMacro(Translations,PointType);
itkGetMacro(Translations,PointType);
itkSetMacro(HiddenTranslations,PointType);
itkGetMacro(HiddenTranslations,PointType);
itkSetMacro(Rotations,PointType);
itkGetMacro(Rotations,PointType);
itkSetMacro(HiddenRotations,PointType);
itkGetMacro(HiddenRotations,PointType);
itkSetMacro(UserTranslations,PointType);
itkGetMacro(UserTranslations,PointType);
@ -484,25 +687,17 @@ public:
itkSetMacro(UserRotations,PointType);
itkGetMacro(UserRotations,PointType);
itkSetMacro(ReferenceTransform,TransformMatrixType);
itkGetMacro(ReferenceTransform,TransformMatrixType);
itkSetMacro(CurrentTransform,TransformMatrixType);
itkGetMacro(CurrentTransform,TransformMatrixType);
PointType GetT();
PointType GetR();
protected:
PointType
m_Translations,
m_Rotations,
m_HiddenTranslations,
m_HiddenRotations,
m_UserTranslations,
m_UserRotations;
TransformMatrixType
m_ReferenceTransform,
m_CurrentTransform;
/** Default Constructor **/
TransformMetaInformation ();

31
reg23Topograms/itkDTRrecon/autoreg/CMakeLists.txt Normal file
View File

@ -0,0 +1,31 @@
SET(LIB_NAME autoreg)
#SET(CMAKE_CXX_FLAGS "-std=c++11 -fPIC")
find_package(ITK REQUIRED)
include(${ITK_USE_FILE})
INCLUDE_DIRECTORIES(${ITK_INCLUDE_DIRS})
SET(HDR
itkDRTHelpers.h
itkgTwoImageToOneImageMetric.h
itkgTwoImageToOneImageMetric.hxx
itkMutualInformationTwoImageToOneImageMetric.h
itkMutualInformationTwoImageToOneImageMetric.hxx
itkNormalizedCorrelationTwoImageToOneImageMetric.h
itkNormalizedCorrelationTwoImageToOneImageMetric.hxx
itkTwoProjectionImageRegistrationMethod.h
itkTwoProjectionImageRegistrationMethod.hxx
)
ADD_LIBRARY(${LIB_NAME} ${HDR})
target_link_libraries(
${LIB_NAME}
${ITK_LIBRARIES}
)
target_include_directories (${LIB_NAME}
PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})

93
reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h Normal file
View File

@ -0,0 +1,93 @@
#ifndef ITKDRTHELPERS_H
#define ITKDRTHELPERS_H
#include "itkMath.h"
#include <cstddef>
#include <cstring>
#include <iterator>
//Function to get method and class name for logging purposes.
template <size_t FL, size_t PFL>
const char* computeMethodName(const char (&function)[FL], const char (&prettyFunction)[PFL])
{
using reverse_ptr = std::reverse_iterator<const char*>;
thread_local static char result[PFL];
const char* locFuncName = std::search(prettyFunction, prettyFunction + PFL - 1, function, function + FL - 1);
const char* locClassName = std::find(reverse_ptr(locFuncName), reverse_ptr(prettyFunction), ' ').base();
const char* endFuncName = std::find(locFuncName, prettyFunction + PFL - 1, '(');
result[0] = '\0';
std::strncat(result, locClassName, endFuncName - locClassName + 1);
std::strcat(result, ")");
return result;
}
#define __COMPACT_PRETTY_FUNCTION__ computeMethodName(__FUNCTION__, __PRETTY_FUNCTION__)
namespace itk {
/* reference string required for comparison with tag values */
static const char* ImageOrientationStrings[] = {
"NotDefined",
"HFS",
"FFS",
"HFP",
"FFP",
};
// constant for converting degrees into radians
const float dtr = itk::Math::pi_over_180;
static const bool verbose = false;
/* this is in the end a IEC to HFS...
* but we keep this for ourselves...
*/
static double Standard_DRT2LPS[9] = {
1, 0, 0,
0, 0, -1,
0, 1, 0
};
static double PAElementsIEC[9] = {
1, 0, 0,
0, -1, 0,
0, 0, -1
};
static double LATElementsIEC[9] = {
0, 0, -1,
0, -1, 0,
-1, 0, 0
};
static double HFS2IEC[9] = {
1, 0, 0,
0, 0, 1,
0, -1, 0
};
static double FFS2IEC[9] = {
-1, 0, 0,
0, 0, -1,
0, -1, 0
};
static double HFP2IEC[9] = {
-1, 0, 0,
0, 0, 1,
0, 1, 0
};
static double FFP2IEC[9] = {
1, 0, 0,
0, 0, -1,
0, 1, 0
};
static double PAT2IEC[9] = {
1, 0, 0,
0, 0, 1,
0, -1, 0
};
}
#endif // ITKDRTHELPERS_H

117
reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.h Normal file
View File

@ -0,0 +1,117 @@
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkMutualInformationTwoImageToOneImageMetric_h
#define itkMutualInformationTwoImageToOneImageMetric_h
#include "itkCovariantVector.h"
#include "itkMacro.h"
#include "itkPoint.h"
#include "itkgTwoImageToOneImageMetric.h"
namespace itk {
/** \class NormalizedCorrelationTwoImageToOneImageMetric
* \brief Computes similarity between two fixed images and one moving image
*
* This metric computes the correlation between pixels in the two fixed images
* and pixels in the moving image. The spatial correspondance between
* two fixed images and the moving image is established through a Transform. Pixel values are
* taken from the fixed images, their positions are mapped to the moving
* image and result in general in non-grid position on it. Values at these
* non-grid position of the moving image are interpolated using user-selected
* Interpolators. The correlation is normalized by the autocorrelations of both
* the fixed and moving images.
*
* \ingroup RegistrationMetrics
* \ingroup TwoProjectionRegistration
*/
template <typename TFixedImage, typename TMovingImage>
class MutualInformationTwoImageToOneImageMetric : public gTwoImageToOneImageMetric<TFixedImage, TMovingImage> {
public:
ITK_DISALLOW_COPY_AND_ASSIGN(MutualInformationTwoImageToOneImageMetric);
/** Standard class type alias. */
using Self = MutualInformationTwoImageToOneImageMetric;
using Superclass = gTwoImageToOneImageMetric<TFixedImage, TMovingImage>;
using Pointer = SmartPointer<Self>;
using ConstPointer = SmartPointer<const Self>;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(MutualInformationTwoImageToOneImageMetric, Object);
/** Types transferred from the base class */
using RealType = typename Superclass::RealType;
using TransformType = typename Superclass::TransformType;
using TransformPointer = typename Superclass::TransformPointer;
using TransformParametersType = typename Superclass::TransformParametersType;
using TransformJacobianType = typename Superclass::TransformJacobianType;
using GradientPixelType = typename Superclass::GradientPixelType;
using MeasureType = typename Superclass::MeasureType;
using DerivativeType = typename Superclass::DerivativeType;
using FixedImageType = typename Superclass::FixedImageType;
using MovingImageType = typename Superclass::MovingImageType;
using FixedImageConstPointer = typename Superclass::FixedImageConstPointer;
using MovingImageConstPointer = typename Superclass::MovingImageConstPointer;
/** Get the derivatives of the match measure. */
void
GetDerivative(const TransformParametersType& parameters, DerivativeType& Derivative) const override;
/** Get the value for single valued optimizers. */
MeasureType
GetValue(const TransformParametersType& parameters) const override;
/** Get the value using the current transforms. */
MeasureType
GetValue() const;
/** Get value and derivatives for multiple valued optimizers. */
void
GetValueAndDerivative(const TransformParametersType& parameters,
MeasureType& Value,
DerivativeType& Derivative) const override;
/** Set/Get SubtractMean boolean. If true, the sample mean is subtracted
* from the sample values in the cross-correlation formula and
* typically results in narrower valleys in the cost fucntion.
* Default value is false. */
itkSetMacro(SubtractMean, bool);
itkGetConstReferenceMacro(SubtractMean, bool);
itkBooleanMacro(SubtractMean);
protected:
MutualInformationTwoImageToOneImageMetric();
~MutualInformationTwoImageToOneImageMetric() override = default;
void
PrintSelf(std::ostream& os, Indent indent) const override;
private:
bool m_SubtractMean;
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkMutualInformationTwoImageToOneImageMetric.hxx"
#endif
#endif

237
reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx Normal file
View File

@ -0,0 +1,237 @@
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkMutualInformationTwoImageToOneImageMetric_hxx
#define itkMutualInformationTwoImageToOneImageMetric_hxx
#include "itkImageRegionConstIteratorWithIndex.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkMutualInformationTwoImageToOneImageMetric.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkNormalizeImageFilter.h"
#include "itkTranslationTransform.h"
#include "itkProgressReporter.h"
#include <limits>
namespace itk {
template <typename TFixedImage, typename TMovingImage>
MutualInformationTwoImageToOneImageMetric<TFixedImage,
TMovingImage>::MutualInformationTwoImageToOneImageMetric()
{
m_SubtractMean = false;
}
template <typename TFixedImage, typename TMovingImage>
typename MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue(
const TransformParametersType& parameters) const
{
FixedImageConstPointer fixedImage1 = this->m_FixedImage1;
if (!fixedImage1) {
itkExceptionMacro(<< "Fixed image1 has not been assigned");
}
FixedImageConstPointer fixedImage2 = this->m_FixedImage2;
if (!fixedImage2) {
itkExceptionMacro(<< "Fixed image2 has not been assigned");
}
if (!this->m_Filter1) {
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
}
if (!this->m_Filter2) {
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
}
this->SetTransformParameters(parameters);
return GetValue();
}
template <typename TFixedImage, typename TMovingImage>
typename MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue() const
{
FixedImageConstPointer fixedImage1 = this->m_FixedImage1;
if (!fixedImage1) {
itkExceptionMacro(<< "Fixed image1 has not been assigned");
}
FixedImageConstPointer fixedImage2 = this->m_FixedImage2;
if (!fixedImage2) {
itkExceptionMacro(<< "Fixed image2 has not been assigned");
}
if (!this->m_Filter1) {
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
}
if (!this->m_Filter2) {
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
}
if (!this->m_Transform1) {
itkExceptionMacro(<< "m_Transform1 has not been assigned");
}
if (!this->m_Transform2) {
itkExceptionMacro(<< "m_Transform2 has not been assigned");
}
// std::cout<< "----- itkMutualInformationTwoImage GetValue -----" <<std::endl;
// std::cout<< "Mutual GetValue T1 :: center : " <<
// this->m_Transform1.GetPointer()->GetCenter() <<std::endl;
// std::cout<< "Mutual GetValue T1 :: Pars : " <<
// this->m_Transform1.GetPointer()->GetParameters() <<std::endl;
// std::cout<< "--------------------------------" <<std::endl;
this->m_Interpolator1->SetTransform(this->m_Transform1);
this->m_Interpolator1->Initialize();
this->m_Filter1->Update();
using InternalMetricType = itk::MattesMutualInformationImageToImageMetric<TFixedImage, TMovingImage>;
/*
* NearestNeighborType in our images results in not resolving for the cranio-caudal direction.
* in LPS-Z we are very much bound to the CT resolution, therefore if we take the Nearest Neighbor
* we simply do not explore that direction... and will always have an uncertainty of 1 slice thickness.
* */
//using NearestNeighborType = itk::NearestNeighborInterpolateImageFunction<TMovingImage, double>;
/*
* LinearInterpolateImageFunction performs better for us, performance degradation not noticeable.
* */
using LinearInterpType = itk::LinearInterpolateImageFunction<TMovingImage, double>;
//We need to set transform and parameters.
auto dummyTransform = TransformType::New();
auto dummyParameters = dummyTransform->GetParameters();
auto dummyInterpolator1 = LinearInterpType::New();
auto metric1 = InternalMetricType::New();
metric1->SetTransform(dummyTransform);
metric1->SetTransformParameters(dummyParameters);
metric1->SetInterpolator(dummyInterpolator1);
//auto fixedImageRegion1 = fixedImage1->GetBufferedRegion();
//metric1->SetFixedImageRegion(fixedImageRegion1);
/* We get the fixed image region from the parent class... */
metric1->SetFixedImageRegion(Superclass::GetFixedImageRegion1());
//std::cout<< "-----> Mutual :: fixedImageRegion1: "<< metric1->GetFixedImageRegion() << std::endl;
auto movingImageRegion = this->m_Filter1->GetOutput()->GetBufferedRegion();
unsigned int numberOfPixels = movingImageRegion.GetNumberOfPixels();
//auto numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
// since we have a ROI, then we should not set allPixels to TRUE.
//metric1->UseAllPixelsOn();
metric1->SetNumberOfHistogramBins(50);
// InternalImageType::IndexType pIndex;
// pIndex[0]=200;
// pIndex[1]=300;
// InternalImageType::PixelType pPixel = fixedImage1->GetPixel(pIndex);
// std::cout<<"MATTES PIXEL: "<<pPixel <<std::endl;
metric1->SetFixedImage(fixedImage1);
metric1->SetMovingImage(this->m_Filter1->GetOutput());
metric1->Initialize();
auto measure1 = metric1->GetValue(dummyParameters);
this->m_Interpolator2->SetTransform(this->m_Transform2);
this->m_Interpolator2->Initialize();
this->m_Filter2->Update();
auto dummyInterpolator2 = LinearInterpType::New();
auto metric2 = InternalMetricType::New();
metric2->SetTransform(dummyTransform);
metric2->SetTransformParameters(dummyParameters);
metric2->SetInterpolator(dummyInterpolator2);
// auto fixedImageRegion2 = fixedImage2->GetBufferedRegion();
// metric2->SetFixedImageRegion(fixedImageRegion2);
metric2->SetFixedImageRegion(Superclass::GetFixedImageRegion2());
//std::cout<< "-----> Mutual :: fixedImageRegion2: "<< metric2->GetFixedImageRegion() << std::endl;
movingImageRegion = this->m_Filter2->GetOutput()->GetBufferedRegion();
numberOfPixels = movingImageRegion.GetNumberOfPixels();
//numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
//metric2->SetNumberOfSpatialSamples(numberOfSamples);
//metric2->UseAllPixelsOn();
metric2->SetNumberOfHistogramBins(50);
metric2->SetFixedImage(fixedImage2);
metric2->SetMovingImage(this->m_Filter2->GetOutput());
metric2->Initialize();
auto measure2 = metric2->GetValue(dummyParameters);
//std::cout << "movingValueMin: " << movingValueMin << " movingValueMax: " << movingValueMax << std::endl;
//std::cout << "fixedValueMin: " << fixedValueMin << " fixedValueMax: " << fixedValueMax << std::endl;
// Calculate the measure value between fixed image 2 and the moving image
auto measure = (measure1 + measure2) / 2.0;
// auto oldprecision = std::cout.precision();
// std::cout.precision(std::numeric_limits<double>::digits10 + 2);
// std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
// std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
// std::cout << "Measure: " << measure << std::endl;
// std::cout << "=======================================================" << std::endl;
// std::cout.precision(oldprecision);
return measure;
}
template <typename TFixedImage, typename TMovingImage>
void MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetDerivative(
const TransformParametersType& itkNotUsed(parameters),
DerivativeType& itkNotUsed(derivative)) const
{
// under construction
}
template <typename TFixedImage, typename TMovingImage>
void MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValueAndDerivative(
const TransformParametersType& itkNotUsed(parameters),
MeasureType& itkNotUsed(value),
DerivativeType& itkNotUsed(derivative)) const
{
// under construction
}
template <typename TFixedImage, typename TMovingImage>
void MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::PrintSelf(std::ostream& os,
Indent indent) const
{
Superclass::PrintSelf(os, indent);
os << indent << "SubtractMean: " << m_SubtractMean << std::endl;
}
} // end namespace itk
#endif

127
reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.h Normal file
View File

@ -0,0 +1,127 @@
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkNormalizedCorrelationTwoImageToOneImageMetric_h
#define itkNormalizedCorrelationTwoImageToOneImageMetric_h
#include "itkCovariantVector.h"
#include "itkMacro.h"
#include "itkNormalizedCorrelationImageToImageMetric.hxx"
#include "itkPoint.h"
#include "itkgTwoImageToOneImageMetric.h"
namespace itk {
/** \class NormalizedCorrelationTwoImageToOneImageMetric
* \brief Computes similarity between two fixed images and one moving image
*
* This metric computes the correlation between pixels in the two fixed images
* and pixels in the moving image. The spatial correspondance between
* two fixed images and the moving image is established through a Transform. Pixel values are
* taken from the fixed images, their positions are mapped to the moving
* image and result in general in non-grid position on it. Values at these
* non-grid position of the moving image are interpolated using user-selected
* Interpolators. The correlation is normalized by the autocorrelations of both
* the fixed and moving images.
*
* \ingroup RegistrationMetrics
* \ingroup TwoProjectionRegistration
*/
template <typename TFixedImage, typename TMovingImage>
class NormalizedCorrelationTwoImageToOneImageMetric : public gTwoImageToOneImageMetric<TFixedImage, TMovingImage> {
public:
ITK_DISALLOW_COPY_AND_ASSIGN(NormalizedCorrelationTwoImageToOneImageMetric);
/** Standard class type alias. */
using Self = NormalizedCorrelationTwoImageToOneImageMetric;
using Superclass = gTwoImageToOneImageMetric<TFixedImage, TMovingImage>;
using Pointer = SmartPointer<Self>;
using ConstPointer = SmartPointer<const Self>;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(NormalizedCorrelationTwoImageToOneImageMetric, Object);
/** Types transferred from the base class */
using RealType = typename Superclass::RealType;
using TransformType = typename Superclass::TransformType;
using TransformPointer = typename Superclass::TransformPointer;
using TransformParametersType = typename Superclass::TransformParametersType;
using TransformJacobianType = typename Superclass::TransformJacobianType;
using GradientPixelType = typename Superclass::GradientPixelType;
using MeasureType = typename Superclass::MeasureType;
using DerivativeType = typename Superclass::DerivativeType;
using FixedImageType = typename Superclass::FixedImageType;
using FixedImageConstPointer = typename Superclass::FixedImageConstPointer;
using FixedImageMaskPointer = typename Superclass::FixedImageMaskPointer;
using FixedImageRegionType = typename Superclass::FixedImageRegionType;
using MovingImageType = typename Superclass::MovingImageType;
using MovingImageConstPointer = typename Superclass::MovingImageConstPointer;
using MovingImageMaskPointer = typename Superclass::MovingImageMaskPointer;
using ChangeInformationFilterPointer = typename Superclass::ChangeInformationFilterPointer;
/** Get the derivatives of the match measure. */
void
GetDerivative(const TransformParametersType& parameters, DerivativeType& Derivative) const override;
/** Get the value for single valued optimizers, after applying paramters to the transforms. */
MeasureType
GetValue(const TransformParametersType& parameters) const override;
/** Get the value using the current transforms. */
MeasureType
GetValue() const;
/** Get value and derivatives for multiple valued optimizers. */
void
GetValueAndDerivative(const TransformParametersType& parameters,
MeasureType& Value,
DerivativeType& Derivative) const override;
/** Set/Get SubtractMean boolean. If true, the sample mean is subtracted
* from the sample values in the cross-correlation formula and
* typically results in narrower valleys in the cost fucntion.
* Default value is false. */
itkSetMacro(SubtractMean, bool);
itkGetConstReferenceMacro(SubtractMean, bool);
itkBooleanMacro(SubtractMean);
protected:
NormalizedCorrelationTwoImageToOneImageMetric();
~NormalizedCorrelationTwoImageToOneImageMetric() override = default;
void
PrintSelf(std::ostream& os, Indent indent) const override;
private:
bool m_SubtractMean;
MeasureType CalculateMeasure(int imageId) const;
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkNormalizedCorrelationTwoImageToOneImageMetric.hxx"
#endif
#endif

317
reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx Normal file
View File

@ -0,0 +1,317 @@
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkNormalizedCorrelationTwoImageToOneImageMetric_hxx
#define itkNormalizedCorrelationTwoImageToOneImageMetric_hxx
#include "itkImageRegionConstIteratorWithIndex.h"
#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h"
#
namespace itk {
template <typename TFixedImage, typename TMovingImage>
NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage,
TMovingImage>::NormalizedCorrelationTwoImageToOneImageMetric()
{
m_SubtractMean = false;
}
template <typename TFixedImage, typename TMovingImage>
typename NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::CalculateMeasure(int imageId) const
{
using FixedIteratorType = itk::ImageRegionConstIteratorWithIndex<FixedImageType>;
using MovingIteratorType = itk::ImageRegionConstIteratorWithIndex<MovingImageType>;
using AccumulateType = typename NumericTraits<MeasureType>::AccumulateType;
FixedImageConstPointer fixedImage;
ChangeInformationFilterPointer filter;
FixedImageMaskPointer fixedImageMask;
FixedImageRegionType fixedImageRegion;
MeasureType measure;
typename FixedImageType::IndexType indexFixed;
typename MovingImageType::IndexType indexMoving;
typename Superclass::InputPointType inputPoint;
if (imageId == 1) {
this->m_Interpolator1->SetTransform(this->m_Transform1);
this->m_Interpolator1->Initialize();
fixedImage = this->m_FixedImage1;
fixedImageMask = this->m_FixedImageMask1;
fixedImageRegion = this->GetFixedImageRegion1();
filter = this->m_Filter1;
} else if (imageId == 2) {
this->m_Interpolator2->SetTransform(this->m_Transform2);
this->m_Interpolator2->Initialize();
fixedImage = this->m_FixedImage2;
fixedImageMask = this->m_FixedImageMask2;
fixedImageRegion = this->GetFixedImageRegion2();
filter = this->m_Filter2;
} else {
itkExceptionMacro(<< "Non supported image id.");
}
FixedIteratorType ti(fixedImage, fixedImageRegion);
this->m_NumberOfPixelsCounted = 0;
int numberOfPixelsVisited = 0;
filter->Update();
MovingImageConstPointer imageDRTIn = filter->GetOutput();
//this->WriteImage(imageDRTIn, "", "");
AccumulateType sff = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType smm = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType sfm = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType sf = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType sm = NumericTraits<AccumulateType>::ZeroValue();
RealType movingValueMin = NumericTraits<RealType>::max();
RealType fixedValueMin = NumericTraits<RealType>::max();
RealType movingValueMax = NumericTraits<RealType>::NonpositiveMin();
RealType fixedValueMax = NumericTraits<RealType>::NonpositiveMin();
RealType maxX = NumericTraits<RealType>::NonpositiveMin();
RealType maxY = NumericTraits<RealType>::NonpositiveMin();
RealType maxZ = NumericTraits<RealType>::NonpositiveMin();
RealType minX = NumericTraits<RealType>::max();
RealType minY = NumericTraits<RealType>::max();
RealType minZ = NumericTraits<RealType>::max();
AccumulateType m_meanCurr_f = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType m_meanCurr_m = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType m_meanPrev_f = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType m_meanPrev_m = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType m_sPrev_f = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType m_sPrev_m = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType m_sCurr_f = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType m_sCurr_m = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType m_varianceCurr_f = NumericTraits<AccumulateType>::ZeroValue();
AccumulateType m_varianceCurr_m = NumericTraits<AccumulateType>::ZeroValue();
while (!ti.IsAtEnd()) {
indexFixed = ti.GetIndex();
++numberOfPixelsVisited;
fixedImage->TransformIndexToPhysicalPoint(indexFixed, inputPoint);
if (fixedImageMask && !fixedImageMask->IsInsideInWorldSpace(inputPoint)) {
++ti;
continue;
}
if (this->m_MovingImageMask && !this->m_MovingImageMask->IsInsideInWorldSpace(inputPoint)) {
++ti;
continue;
}
const RealType fixedValue = ti.Get();
if (fixedValue == -1024) {
++ti;
continue;
}
if (imageDRTIn->TransformPhysicalPointToIndex(inputPoint, indexMoving)) {
double x = inputPoint[0];
double y = inputPoint[1];
double z = inputPoint[2];
maxX = maxX > x ? maxX : x;
minX = minX < x ? minX : x;
maxY = maxY > y ? maxY : y;
minY = minY < y ? minY : y;
maxZ = maxZ > z ? maxZ : z;
minZ = minZ < z ? minZ : z;
const RealType movingValue = imageDRTIn->GetPixel(indexMoving);
if (NumericTraits<short>::NonpositiveMin() == movingValue) {
++ti;
continue;
}
movingValueMin = movingValue < movingValueMin ? movingValue : movingValueMin;
movingValueMax = movingValue > movingValueMax ? movingValue : movingValueMax;
fixedValueMin = fixedValue < fixedValueMin ? fixedValue : fixedValueMin;
fixedValueMax = fixedValue > fixedValueMax ? fixedValue : fixedValueMax;
sff += fixedValue * fixedValue;
smm += movingValue * movingValue;
sfm += fixedValue * movingValue;
if (this->m_SubtractMean) {
sf += fixedValue;
sm += movingValue;
}
this->m_NumberOfPixelsCounted++;
if (this->m_NumberOfPixelsCounted == 1) {
m_meanCurr_f = fixedValue;
m_meanCurr_m = movingValue;
} else {
m_meanPrev_f = m_meanCurr_f;
m_meanPrev_m = m_meanCurr_m;
m_sPrev_f = m_sCurr_f;
m_sPrev_m = m_sCurr_m;
m_meanCurr_f = m_meanPrev_f + (fixedValue - m_meanPrev_f) / this->m_NumberOfPixelsCounted;
m_sCurr_f = m_sPrev_f + (fixedValue - m_meanPrev_f) * (fixedValue - m_meanCurr_f);
m_varianceCurr_f = m_sCurr_f / (this->m_NumberOfPixelsCounted - 1);
m_meanCurr_m = m_meanPrev_m + (movingValue - m_meanPrev_m) / this->m_NumberOfPixelsCounted;
m_sCurr_m = m_sPrev_m + (movingValue - m_meanPrev_m) * (movingValue - m_meanCurr_m);
m_varianceCurr_m = m_sCurr_m / (this->m_NumberOfPixelsCounted - 1);
}
}
++ti;
}
if (this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0) {
sff -= (sf * sf / this->m_NumberOfPixelsCounted);
smm -= (sm * sm / this->m_NumberOfPixelsCounted);
sfm -= (sf * sm / this->m_NumberOfPixelsCounted);
}
RealType denom = -1.0 * sqrt(sff * smm);
if (this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
measure = sfm / denom;
} else {
measure = NumericTraits<MeasureType>::Zero;
}
// std::cout << "movingValueMin: " << movingValueMin << " movingValueMax: " << movingValueMax << std::endl;
// std::cout << "fixedValueMin: " << fixedValueMin << " fixedValueMax: " << fixedValueMax << std::endl;
// std::cout << "max coordinates: " << maxX << " " << maxX << " " << maxZ << std::endl;
// std::cout << "min coordinates: " << minX << " " << minY << " " << minZ << std::endl;
return measure;
}
template <typename TFixedImage, typename TMovingImage>
typename NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue(
const TransformParametersType& parameters) const
{
if (!this->m_FixedImage1) {
itkExceptionMacro(<< "Fixed image1 has not been assigned");
}
if (!this->m_FixedImage2) {
itkExceptionMacro(<< "Fixed image2 has not been assigned");
}
if (!this->m_Filter1) {
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
}
if (!this->m_Filter2) {
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
}
if (!this->SetTransformParameters(parameters)) {
std::cout << "Measure: 0.0, transform Invalid!" << std::endl;
std::cout << "=======================================================" << std::endl;
return 0.0;
}
return GetValue();
}
template <typename TFixedImage, typename TMovingImage>
typename NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::MeasureType
NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue() const
{
if (!this->m_FixedImage1) {
itkExceptionMacro(<< "Fixed image1 has not been assigned");
}
if (!this->m_FixedImage2) {
itkExceptionMacro(<< "Fixed image2 has not been assigned");
}
if (!this->m_Filter1) {
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
}
if (!this->m_Filter2) {
itkExceptionMacro(<< "Fixed filter1 has not been assigned");
}
// Calculate the measure value between fixed image 1 and the moving image
auto oldprecision = std::cout.precision();
// std::cout<<"Region " <<this->GetFixedImageRegion1() <<std::endl;
// std::cout<<"Buffered Region " <<this->GetFixedImage1()->GetBufferedRegion() <<std::endl;
// std::cout<<"Buffered Region " <<this->m_Filter1->GetOutput()->GetBufferedRegion() <<std::endl;
MeasureType measure1 = CalculateMeasure(1);
// std::cout.precision(std::numeric_limits<double>::digits10 + 2);
// std::cout << "Measure1: " << measure1 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
// std::cout.precision(oldprecision);
MeasureType measure2 = CalculateMeasure(2);
// std::cout.precision(std::numeric_limits<double>::digits10 + 2);
// std::cout << "Measure2: " << measure2 << " Counts: " << this->m_NumberOfPixelsCounted << std::endl;
auto measure = (measure1 + measure2) / 2.0;
// std::cout << "Measure: " << measure << std::endl;
// std::cout << "=======================================================" << std::endl;
// std::cout.precision(oldprecision);
return measure;
}
template <typename TFixedImage, typename TMovingImage>
void NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetDerivative(
const TransformParametersType& itkNotUsed(parameters),
DerivativeType& itkNotUsed(derivative)) const
{
// under construction
}
template <typename TFixedImage, typename TMovingImage>
void NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValueAndDerivative(
const TransformParametersType& itkNotUsed(parameters),
MeasureType& itkNotUsed(value),
DerivativeType& itkNotUsed(derivative)) const
{
// under construction
}
template <typename TFixedImage, typename TMovingImage>
void NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::PrintSelf(std::ostream& os,
Indent indent) const
{
Superclass::PrintSelf(os, indent);
os << indent << "SubtractMean: " << m_SubtractMean << std::endl;
}
} // end namespace itk
#endif

294
reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.h Normal file
View File

@ -0,0 +1,294 @@
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkTwoProjectionImageRegistrationMethod_h
#define itkTwoProjectionImageRegistrationMethod_h
#include "itkChangeInformationImageFilter.h"
#include "itkDataObjectDecorator.h"
#include "itkImage.h"
#include "itkMacro.h"
#include "itkProcessObject.h"
#include "itkSingleValuedNonLinearOptimizer.h"
#include "itkgTwoImageToOneImageMetric.h"
#include "itkObject.h"
#include "itkEventObject.h"
#include <iostream>
namespace itk {
/** \class TwoProjectionImageRegistrationMethod
* \brief Base class for Image Registration Methods
*
* This Class define the generic interface for a registration method.
*
* This class is templated over the type of the two image to be
* registered. A generic Transform is used by this class. That allows
* to select at run time the particular type of transformation that
* is to be applied for registering the images.
*
* This method use a generic Metric in order to compare the two images.
* the final goal of the registration method is to find the set of
* parameters of the Transformation that optimizes the metric.
*
* The registration method also support a generic optimizer that can
* be selected at run-time. The only restriction for the optimizer is
* that it should be able to operate in single-valued cost functions
* given that the metrics used to compare images provide a single
* value as output.
*
* The terms : Fixed image and Moving image are used in this class
* to indicate what image is being mapped by the transform.
*
* This class uses the coordinate system of the Fixed image as a reference
* and searchs for a Transform that will map points from the space of the
* Fixed image to the space of the Moving image.
*
* For doing so, a Metric will be continously applied to compare the Fixed
* image with the Transformed Moving image. This process also requires to
* interpolate values from the Moving image.
*
* \ingroup RegistrationFilters
* \ingroup TwoProjectionRegistration
*/
template <typename TFixedImage, typename TMovingImage>
class TwoProjectionImageRegistrationMethod : public ProcessObject {
public:
ITK_DISALLOW_COPY_AND_ASSIGN(TwoProjectionImageRegistrationMethod);
/** Standard class type alias. */
using Self = TwoProjectionImageRegistrationMethod;
using Superclass = ProcessObject;
using Pointer = SmartPointer<Self>;
using ConstPointer = SmartPointer<const Self>;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(TwoProjectionImageRegistrationMethod, ProcessObject);
/** Type of the Fixed image. */
using FixedImageType = TFixedImage;
using FixedImageConstPointer = typename FixedImageType::ConstPointer;
/** Type of the Moving image. */
using MovingImageType = TMovingImage;
using MovingImageConstPointer = typename MovingImageType::ConstPointer;
using ChangeInformationFilterType = ChangeInformationImageFilter<MovingImageType>;
using ChangeInformationFilterPointer = typename ChangeInformationFilterType::Pointer;
/** Type of the metric. */
using MetricType = gTwoImageToOneImageMetric<FixedImageType, MovingImageType>;
using MetricPointer = typename MetricType::Pointer;
using FixedImageRegionType = typename MetricType::FixedImageRegionType;
/** Type of the Transform . */
using TransformType = typename MetricType::TransformType;
using TransformPointer = typename TransformType::Pointer;
/** Type for the output: Using Decorator pattern for enabling
* the Transform to be passed in the data pipeline */
using TransformOutputType = DataObjectDecorator<TransformType>;
using TransformOutputPointer = typename TransformOutputType::Pointer;
using TransformOutputConstPointer = typename TransformOutputType::ConstPointer;
/** Type of the Interpolator. */
using InterpolatorType = typename MetricType::InterpolatorType;
using InterpolatorPointer = typename InterpolatorType::Pointer;
/** Type of the optimizer. */
using OptimizerType = SingleValuedNonLinearOptimizer;
/** Type of the Transformation parameters This is the same type used to
* represent the search space of the optimization algorithm */
using ParametersType = typename MetricType::TransformParametersType;
/** Smart Pointer type to a DataObject. */
using DataObjectPointer = typename DataObject::Pointer;
/** Method that initiates the registration. This will Initialize and ensure
* that all inputs the registration needs are in place, via a call to
* Initialize() will then start the optimization process via a call to
* StartOptimization() */
void
StartRegistration();
/** Method that initiates the optimization process. */
void
StartOptimization();
/** Set/Get the Fixed images. */
void
SetFixedImage1(const FixedImageType* fixedImage1);
void
SetFixedImage2(const FixedImageType* fixedImage2);
itkGetConstObjectMacro(FixedImage1, FixedImageType);
itkGetConstObjectMacro(FixedImage2, FixedImageType);
/** Set/Get the Moving image. */
void
SetMovingImage(const MovingImageType* movingImage);
itkGetConstObjectMacro(MovingImage, MovingImageType);
/** Set/Get the Optimizer. */
itkSetObjectMacro(Optimizer, OptimizerType);
itkGetConstObjectMacro(Optimizer, OptimizerType);
/** Set/Get the Metric. */
itkSetObjectMacro(Metric, MetricType);
itkGetConstObjectMacro(Metric, MetricType);
/** Set/Get the Transfrom1. */
itkSetObjectMacro(Transform1, TransformType);
itkGetConstObjectMacro(Transform1, TransformType);
/** Set/Get the Transfrom2. */
itkSetObjectMacro(Transform2, TransformType);
itkGetConstObjectMacro(Transform2, TransformType);
/** Set/Get the IsocIECTransform. */
itkSetObjectMacro(IsocIECTransform, TransformType);
itkGetConstObjectMacro(IsocIECTransform, TransformType);
/** Set/Get the Interpolators. */
itkSetObjectMacro(Interpolator1, InterpolatorType);
itkSetObjectMacro(Interpolator2, InterpolatorType);
itkGetConstObjectMacro(Interpolator1, InterpolatorType);
itkGetConstObjectMacro(Interpolator2, InterpolatorType);
/** Set/Get the Meta informations. */
itkSetObjectMacro(TransformMetaInfo, R23MetaInformation);
itkGetConstObjectMacro(TransformMetaInfo, R23MetaInformation);
itkSetObjectMacro(internalMeta1, InternalTransformMetaInformation);
itkGetConstObjectMacro(internalMeta1, InternalTransformMetaInformation);
itkSetObjectMacro(internalMeta2, InternalTransformMetaInformation);
itkGetConstObjectMacro(internalMeta2, InternalTransformMetaInformation);
/** Set/Get the output filters. */
itkSetObjectMacro(Filter1, ChangeInformationFilterType);
itkGetConstObjectMacro(Filter1, ChangeInformationFilterType);
itkSetObjectMacro(Filter2, ChangeInformationFilterType);
itkGetConstObjectMacro(Filter2, ChangeInformationFilterType);
/** Set the region of the fixed image to be considered as region of
interest during the registration. This region will be passed to
the ImageMetric in order to restrict the metric computation to
consider only this region.
\warning The same region can also be set directly into the metric.
please avoid to set the region in both places since this can lead
to inconsistent configurations. */
void
SetFixedImageRegion1(const FixedImageRegionType& region1);
void
SetFixedImageRegion2(const FixedImageRegionType& region2);
/** Get the region of the fixed image to be considered as region of
interest during the registration. This region will be passed to
the ImageMetric in order to restrict the metric computation to
consider only this region. */
itkGetConstReferenceMacro(FixedImageRegion1, FixedImageRegionType);
itkGetConstReferenceMacro(FixedImageRegion2, FixedImageRegionType);
/** True if a region has been defined for the fixed image to which
the ImageMetric will limit its computation */
itkGetMacro(FixedImageRegionDefined1, bool);
itkGetMacro(FixedImageRegionDefined2, bool);
/** Turn on/off the use of a fixed image region to which
the ImageMetric will limit its computation.
\warning The region must have been previously defined using the
SetFixedImageRegion member function */
itkSetMacro(FixedImageRegionDefined1, bool);
itkSetMacro(FixedImageRegionDefined2, bool);
/** Initialize by setting the interconnects between the components. */
virtual void
Initialize();
/** Returns the transform resulting from the registration process */
const TransformOutputType*
GetOutput() const;
/** Make a DataObject of the correct type to be used as the specified
* output. */
using Superclass::MakeOutput;
DataObjectPointer
MakeOutput(DataObjectPointerArraySizeType idx) override;
protected:
TwoProjectionImageRegistrationMethod();
~TwoProjectionImageRegistrationMethod() override = default;
void
PrintSelf(std::ostream& os, Indent indent) const override;
/** Method invoked by the pipeline in order to trigger the computation of
* the registration. */
void
GenerateData() override;
/** Provides derived classes with the ability to set this private var */
itkSetMacro(LastOptimizerParameters, ParametersType);
/** Entry-point for internal ITK filter observer. **/
// void OnInternalFilterProgressReceptor(itk::Object *caller,
// const itk::EventObject &event);
private:
MetricPointer m_Metric;
OptimizerType::Pointer m_Optimizer;
MovingImageConstPointer m_MovingImage;
FixedImageConstPointer m_FixedImage1;
FixedImageConstPointer m_FixedImage2;
TransformPointer m_IsocIECTransform;
TransformPointer m_Transform1;
TransformPointer m_Transform2;
InterpolatorPointer m_Interpolator1;
InterpolatorPointer m_Interpolator2;
InternalTransformMetaInformation::Pointer
m_internalMeta1,
m_internalMeta2;
R23MetaInformation::Pointer
m_TransformMetaInfo;
ChangeInformationFilterPointer
m_Filter1,
m_Filter2;
ParametersType m_InitialOptimizerParameters;
ParametersType m_LastOptimizerParameters;
bool m_FixedImageRegionDefined1;
bool m_FixedImageRegionDefined2;
FixedImageRegionType m_FixedImageRegion1;
FixedImageRegionType m_FixedImageRegion2;
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkTwoProjectionImageRegistrationMethod.hxx"
#endif
#endif

381
reg23Topograms/itkDTRrecon/autoreg/itkTwoProjectionImageRegistrationMethod.hxx Normal file
View File

@ -0,0 +1,381 @@
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkTwoProjectionImageRegistrationMethod_hxx
#define itkTwoProjectionImageRegistrationMethod_hxx
#include "itkCommand.h"
#include "itkTwoProjectionImageRegistrationMethod.h"
namespace itk {
// TwoProjectionImageRegistrationMethod from Jian Wu.
// This part is not really needed and functionality should be moved to
// to the itkImageProcessor.cpp. This is just abstraction with an additional step.
template <typename TFixedImage, typename TMovingImage>
TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::TwoProjectionImageRegistrationMethod()
{
this->SetNumberOfRequiredOutputs(1); // for the Transform
m_FixedImage1 = nullptr; // has to be provided by the user.
m_FixedImage2 = nullptr; // has to be provided by the user.
m_MovingImage = nullptr; // has to be provided by the user.
m_Transform1 = nullptr; // has to be provided by the user.
m_Transform2 = nullptr; // has to be provided by the user.
m_IsocIECTransform = nullptr;
m_TransformMetaInfo = nullptr; // has to be provided by the user.
m_internalMeta1 = nullptr;
m_internalMeta2 = nullptr;
m_Interpolator1 = nullptr; // has to be provided by the user.
m_Interpolator2 = nullptr; // has to be provided by the user.
m_Metric = nullptr; // has to be provided by the user.
m_Optimizer = nullptr; // has to be provided by the user.
m_InitialOptimizerParameters = ParametersType(1);
m_LastOptimizerParameters = ParametersType(1);
m_InitialOptimizerParameters.Fill(0.0f);
m_LastOptimizerParameters.Fill(0.0f);
m_FixedImageRegionDefined1 = false;
m_FixedImageRegionDefined2 = false;
TransformOutputPointer transformDecorator = static_cast<TransformOutputType*>(this->MakeOutput(0).GetPointer());
this->ProcessObject::SetNthOutput(0, transformDecorator.GetPointer());
}
/*
* Set the region of the fixed image 1 to be considered for registration
*/
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedImageRegion1(
const FixedImageRegionType& region1)
{
m_FixedImageRegion1 = region1;
m_FixedImageRegionDefined1 = true;
}
/*
* Set the region of the fixed image 2 to be considered for registration
*/
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedImageRegion2(
const FixedImageRegionType& region2)
{
m_FixedImageRegion2 = region2;
m_FixedImageRegionDefined2 = true;
}
/*
* Initialize by setting the interconnects between components.
*/
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::Initialize()
{
if (!m_FixedImage1) {
itkExceptionMacro(<< "FixedImage1 is not present");
}
if (!m_FixedImage2) {
itkExceptionMacro(<< "FixedImage2 is not present");
}
if (!m_MovingImage) {
itkExceptionMacro(<< "MovingImage is not present");
}
if (!m_Metric) {
itkExceptionMacro(<< "Metric is not present");
}
if (!m_Optimizer) {
itkExceptionMacro(<< "Optimizer is not present");
}
if (!m_Transform1) {
itkExceptionMacro(<< "Transform1 is not present");
}
if (!m_Transform2) {
itkExceptionMacro(<< "Transform2 is not present");
}
//
// Connect the transform to the Decorator.
auto* transformOutput = static_cast<TransformOutputType*>(this->ProcessObject::GetOutput(0));
//transformOutput->Set(m_Transform1.GetPointer());
transformOutput->Set(m_IsocIECTransform.GetPointer());
if (!m_Interpolator1) {
itkExceptionMacro(<< "Interpolator1 is not present");
}
if (!m_Interpolator2) {
itkExceptionMacro(<< "Interpolator2 is not present");
}
if (!m_Filter1) {
itkExceptionMacro(<< "Filter1 is not present");
}
if (!m_Filter1) {
itkExceptionMacro(<< "Filter2 is not present");
}
if (!m_TransformMetaInfo) {
itkExceptionMacro(<< "TransformMetaInfo is not present");
}
if (!m_internalMeta1) {
itkExceptionMacro(<< "internalMeta1 is not present");
}
if (!m_internalMeta2) {
itkExceptionMacro(<< "internalMeta1 is not present");
}
// Setup the metric
m_Metric->SetMovingImage(m_MovingImage);
m_Metric->SetFixedImage1(m_FixedImage1);
m_Metric->SetFixedImage2(m_FixedImage2);
m_Metric->SetTransform1(m_Transform1);
m_Metric->SetTransform2(m_Transform2);
m_Metric->SetIsocTransf(m_IsocIECTransform);
m_Metric->SetTransformMetaInfo(m_TransformMetaInfo);
m_Metric->SetinternalMeta1(m_internalMeta1);
m_Metric->SetinternalMeta2(m_internalMeta2);
m_Metric->SetFilter1(m_Filter1);
m_Metric->SetFilter2(m_Filter2);
m_Metric->SetInterpolator1(m_Interpolator1);
m_Metric->SetInterpolator2(m_Interpolator2);
if (m_FixedImageRegionDefined1) {
m_Metric->SetFixedImageRegion1(m_FixedImageRegion1);
std::cout << "m_FixedImageRegionDefined1 is defined " << std::endl;
} else {
m_Metric->SetFixedImageRegion1(m_FixedImage1->GetBufferedRegion());
}
if (m_FixedImageRegionDefined2) {
m_Metric->SetFixedImageRegion2(m_FixedImageRegion2);
std::cout << "m_FixedImageRegionDefined2 is defined " << std::endl;
} else {
m_Metric->SetFixedImageRegion2(m_FixedImage2->GetBufferedRegion());
}
m_Metric->Initialize();
// Setup the optimizer
m_Optimizer->SetScales(m_Metric->GetWeightings());
m_Optimizer->SetCostFunction(m_Metric);
m_InitialOptimizerParameters =
m_Metric->GetParameters();
m_Optimizer->SetInitialPosition(m_InitialOptimizerParameters);
}
//template <typename TFixedImage, typename TMovingImage>
//void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::OnInternalFilterProgressReceptor(
// itk::Object *caller, const itk::EventObject &event)
//{
// itk::ProcessObject *filter = dynamic_cast<itk::ProcessObject *>(caller);
// if(!itk::ProgressEvent().CheckEvent(&event) || !filter)
// return;
// if (filter)
// {
//// double p = m_CurrentProgressStart;
//// // filter->GetProgress() in [0;1]
//// p += m_CurrentProgressSpan * filter->GetProgress();
//// // NOTE: filter is buggy - not in [0;1] if multi-threading is activated!
//// if (p > (m_CurrentProgressStart + m_CurrentProgressSpan))
//// p = m_CurrentProgressStart + m_CurrentProgressSpan;
//// TaskProgressInfo(m_CurrentProgressDirection, p);
//// if (m_CancelRequest)
//// filter->SetAbortGenerateData(true); // may be handled by filter
// std::cout << "OnInternalFilterProgressReceptor :: filter (TRUE) " << std::endl;
// }
//}
/*
* Starts the Registration Process
*/
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::StartRegistration()
{
itk::ProgressReporter progress(
this,
1, 1000,100);
ParametersType empty(1);
empty.Fill(0.0);
try {
// initialize the interconnects between components
this->Initialize();
if (GetDebug()) {
m_Metric->Print(std::cout);
m_Optimizer->Print(std::cout);
this->Print(std::cout);
}
} catch (ExceptionObject& err) {
m_LastOptimizerParameters = empty;
// pass exception to caller
throw err;
}
this->StartOptimization();
}
/*
* Starts the Optimization process
*/
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::StartOptimization()
{
try {
// do the optimization
m_Optimizer->StartOptimization();
} catch (ExceptionObject& err) {
// An error has occurred in the optimization.
// Update the parameters
m_LastOptimizerParameters = m_Optimizer->GetCurrentPosition();
// Pass exception to caller
throw err;
}
// get the results
m_LastOptimizerParameters = m_Optimizer->GetCurrentPosition();
//Update transform1 and transform2 with the last transform parameters from the optimizer.
m_Metric->SetTransformParameters(m_LastOptimizerParameters);
}
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::PrintSelf(std::ostream& os, Indent indent) const
{
Superclass::PrintSelf(os, indent);
os << indent << "Metric: " << m_Metric.GetPointer() << std::endl;
os << indent << "Optimizer: " << m_Optimizer.GetPointer() << std::endl;
os << indent << "Transform1: " << m_Transform1.GetPointer() << std::endl;
os << indent << "Transform2: " << m_Transform2.GetPointer() << std::endl;
os << indent << "Interpolator 1: " << m_Interpolator1.GetPointer() << std::endl;
os << indent << "Interpolator 2: " << m_Interpolator2.GetPointer() << std::endl;
os << indent << "Fixed Image 1: " << m_FixedImage1.GetPointer() << std::endl;
os << indent << "Fixed Image 2: " << m_FixedImage2.GetPointer() << std::endl;
os << indent << "Moving Image: " << m_MovingImage.GetPointer() << std::endl;
os << indent << "Transform MetaInfo: " << m_TransformMetaInfo.GetPointer() << std::endl;
os << indent << "Fixed Image 1 Region Defined: " << m_FixedImageRegionDefined1 << std::endl;
os << indent << "Fixed Image 2 Region Defined: " << m_FixedImageRegionDefined2 << std::endl;
os << indent << "Fixed Image 1 Region: " << m_FixedImageRegion1 << std::endl;
os << indent << "Fixed Image 2 Region: " << m_FixedImageRegion2 << std::endl;
os << indent << "Initial Optimizer Parameters: " << m_InitialOptimizerParameters << std::endl;
os << indent << "Last Optimizer Parameters: " << m_LastOptimizerParameters << std::endl;
}
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::GenerateData()
{
this->StartRegistration();
}
template <typename TFixedImage, typename TMovingImage>
const typename TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::TransformOutputType*
TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::GetOutput() const
{
return static_cast<const TransformOutputType*>(this->ProcessObject::GetOutput(0));
}
template <typename TFixedImage, typename TMovingImage>
DataObject::Pointer
TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::MakeOutput(DataObjectPointerArraySizeType output)
{
switch (output) {
case 0:
return static_cast<DataObject*>(TransformOutputType::New().GetPointer());
break;
default:
itkExceptionMacro("MakeOutput request for an output number larger than the expected number of outputs");
return nullptr;
}
}
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedImage1(const FixedImageType* fixedImage1)
{
itkDebugMacro("setting Fixed Image 1 to " << fixedImage1);
if (this->m_FixedImage1.GetPointer() != fixedImage1) {
this->m_FixedImage1 = fixedImage1;
// Process object is not const-correct so the const_cast is required here
this->ProcessObject::SetNthInput(0, const_cast<FixedImageType*>(fixedImage1));
this->Modified();
}
}
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetFixedImage2(const FixedImageType* fixedImage2)
{
itkDebugMacro("setting Fixed Image 2 to " << fixedImage2);
if (this->m_FixedImage2.GetPointer() != fixedImage2) {
this->m_FixedImage2 = fixedImage2;
// Process object is not const-correct so the const_cast is required here
this->ProcessObject::SetNthInput(0, const_cast<FixedImageType*>(fixedImage2));
this->Modified();
}
}
template <typename TFixedImage, typename TMovingImage>
void TwoProjectionImageRegistrationMethod<TFixedImage, TMovingImage>::SetMovingImage(const MovingImageType* movingImage)
{
itkDebugMacro("setting Moving Image to " << movingImage);
if (this->m_MovingImage.GetPointer() != movingImage) {
this->m_MovingImage = movingImage;
// Process object is not const-correct so the const_cast is required here
this->ProcessObject::SetNthInput(1, const_cast<MovingImageType*>(movingImage));
this->Modified();
}
}
} // end namespace itk
#endif

319
reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.h Normal file
View File

@ -0,0 +1,319 @@
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkgTwoImageToOneImageMetric_h
#define itkgTwoImageToOneImageMetric_h
#include "DRTMetaInformation.h"
#include "itkChangeInformationImageFilter.h"
#include "itkDRTHelpers.h"
#include "itkGradientRecursiveGaussianImageFilter.h"
#include "itkImageBase.h"
#include "itkInterpolateImageFunction.h"
#include "itkOptimizer.h"
#include "itkResampleImageFilter.h"
#include "itkSingleValuedCostFunction.h"
#include "itkSpatialObject.h"
#include "itkTransform.h"
#include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
#include "itkImageFileWriter.h"
#include "itkRescaleIntensityImageFilter.h"
#include <cstring>
#include <iostream>
#include <libgen.h>
#include <linux/limits.h>
#include <unistd.h>
#include <itkImageProcessorHelpers.h>
namespace itk {
/** \class gTwoImageToOneImageMetric
* \brief Computes similarity between two fixed images and one fixed image.
*
* This Class is templated over the type of the two input images.
* It expects a Transform and two Interpolators to be plugged in.
* This particular class is the base class for a hierarchy of
* similarity metrics.
*
* This class computes a value that measures the similarity
* between two Fixed image and the transformed Moving images.
* The Interpolators are used to compute intensity values on
* non-grid positions resulting from mapping points through
* the Transform.
*
*
* \ingroup RegistrationMetrics
* \ingroup TwoProjectionRegistration
*
*/
template <typename TFixedImage, typename TMovingImage>
class gTwoImageToOneImageMetric : public SingleValuedCostFunction {
public:
ITK_DISALLOW_COPY_AND_ASSIGN(gTwoImageToOneImageMetric);
/** Standard class type alias. */
using Self = gTwoImageToOneImageMetric;
using Superclass = SingleValuedCostFunction;
using Pointer = SmartPointer<Self>;
using ConstPointer = SmartPointer<const Self>;
/** Type used for representing point components */
using CoordinateRepresentationType = Superclass::ParametersValueType;
/** Run-time type information (and related methods). */
itkTypeMacro(gTwoImageToOneImageMetric, SingleValuedCostFunction);
/** Type of the moving Image. */
using MovingImageType = TMovingImage;
using MovingImagePixelType = typename TMovingImage::PixelType;
using MovingImageConstPointer = typename MovingImageType::ConstPointer;
/** Type of the fixed Image. */
using FixedImageType = TFixedImage;
using FixedImageConstPointer = typename FixedImageType::ConstPointer;
using FixedImageRegionType = typename FixedImageType::RegionType;
/** Constants for the image dimensions */
static constexpr unsigned int MovingImageDimension = TMovingImage::ImageDimension;
static constexpr unsigned int FixedImageDimension = TFixedImage::ImageDimension;
/** Type of the Transform Base class */
using TransformType = itk::Euler3DTransform<double>;
using TransformPointer = typename TransformType::Pointer;
using InputPointType = typename TransformType::InputPointType;
using OutputPointType = typename TransformType::OutputPointType;
using TransformParametersType = typename TransformType::ParametersType;
using TransformJacobianType = typename TransformType::JacobianType;
/** Type of the Interpolator Base class */
using InterpolatorType = gSiddonJacobsRayCastInterpolateImageFunction<MovingImageType, CoordinateRepresentationType>;
/** Gaussian filter to compute the gradient of the Moving Image */
using RealType = typename NumericTraits<MovingImagePixelType>::RealType;
using GradientPixelType = CovariantVector<RealType, itkGetStaticConstMacro(MovingImageDimension)>;
using GradientImageType = Image<GradientPixelType, itkGetStaticConstMacro(MovingImageDimension)>;
using GradientImagePointer = SmartPointer<GradientImageType>;
using GradientImageFilterType = GradientRecursiveGaussianImageFilter<MovingImageType, GradientImageType>;
using GradientImageFilterPointer = typename GradientImageFilterType::Pointer;
using InterpolatorPointer = typename InterpolatorType::Pointer;
using ResampleImageFilterType = ResampleImageFilter<MovingImageType, MovingImageType>;
using ResampleImageFilterPointer = typename ResampleImageFilterType::Pointer;
using ChangeInformationFilterType = ChangeInformationImageFilter<MovingImageType>;
using ChangeInformationFilterPointer = typename ChangeInformationFilterType::Pointer;
/** Type for the mask of the fixed image. Only pixels that are "inside"
this mask will be considered for the computation of the metric */
typedef SpatialObject<itkGetStaticConstMacro(FixedImageDimension)> FixedImageMaskType;
using FixedImageMaskPointer = typename FixedImageMaskType::Pointer;
/** Type for the mask of the moving image. Only pixels that are "inside"
this mask will be considered for the computation of the metric */
typedef SpatialObject<itkGetStaticConstMacro(MovingImageDimension)> MovingImageMaskType;
using MovingImageMaskPointer = typename MovingImageMaskType::Pointer;
/** Type of the measure. */
using MeasureType = Superclass::MeasureType;
/** Type of the derivative. */
using DerivativeType = Superclass::DerivativeType;
/** Type of the Transformation parameters This is the same type used to
* represent the search space of the optimization algorithm */
using ParametersType = Superclass::ParametersType;
/** Connect the Fixed Image. */
itkSetConstObjectMacro(FixedImage1, FixedImageType);
/** Connect the Fixed Image. */
itkSetConstObjectMacro(FixedImage2, FixedImageType);
/** Get the Fixed Image. */
itkGetConstObjectMacro(FixedImage1, FixedImageType);
/** Get the Fixed Image. */
itkGetConstObjectMacro(FixedImage2, FixedImageType);
/** Connect the Moving Image. */
itkSetConstObjectMacro(MovingImage, MovingImageType);
/** Get the Moving Image. */
itkGetConstObjectMacro(MovingImage, MovingImageType);
/** Connect the Transform1. */
itkSetObjectMacro(Transform1, TransformType);
/** Get a pointer to the Transform1. */
itkGetConstObjectMacro(Transform1, TransformType);
/** Connect the Transform2. */
itkSetObjectMacro(Transform2, TransformType);
/** Get a pointer to the Transform2. */
itkGetConstObjectMacro(Transform2, TransformType);
itkSetObjectMacro(IsocTransf, TransformType);
itkGetObjectMacro(IsocTransf, TransformType);
/** Connect the Interpolator. */
itkSetObjectMacro(Interpolator1, InterpolatorType);
/** Get a pointer to the Interpolator. */
itkGetConstObjectMacro(Interpolator1, InterpolatorType);
/** Connect the Interpolator. */
itkSetObjectMacro(Interpolator2, InterpolatorType);
/** Get a pointer to the Interpolator. */
itkGetConstObjectMacro(Interpolator2, InterpolatorType);
/** Get the number of pixels considered in the computation. */
itkGetConstReferenceMacro(NumberOfPixelsCounted, unsigned long);
/** Connect the DRTGeometryMetaInfo. */
itkSetObjectMacro(TransformMetaInfo, R23MetaInformation);
/** Get a pointer to the DRTGeometryMetaInfo. */
itkGetConstObjectMacro(TransformMetaInfo, R23MetaInformation);
itkSetObjectMacro(internalMeta1, InternalTransformMetaInformation);
itkGetConstObjectMacro(internalMeta1, InternalTransformMetaInformation);
itkSetObjectMacro(internalMeta2, InternalTransformMetaInformation);
itkGetConstObjectMacro(internalMeta2, InternalTransformMetaInformation);
/** Set the region over which the metric will be computed */
itkSetMacro(FixedImageRegion1, FixedImageRegionType);
/** Set the region over which the metric will be computed */
itkSetMacro(FixedImageRegion2, FixedImageRegionType);
/** Get the region over which the metric will be computed */
itkGetConstReferenceMacro(FixedImageRegion1, FixedImageRegionType);
/** Get the region over which the metric will be computed */
itkGetConstReferenceMacro(FixedImageRegion2, FixedImageRegionType);
/** Set/Get the output filters. */
itkSetObjectMacro(Filter1, ChangeInformationFilterType);
itkGetConstObjectMacro(Filter1, ChangeInformationFilterType);
itkSetObjectMacro(Filter2, ChangeInformationFilterType);
itkGetConstObjectMacro(Filter2, ChangeInformationFilterType);
/** Set/Get the moving image mask. */
itkSetObjectMacro(MovingImageMask, MovingImageMaskType);
itkGetConstObjectMacro(MovingImageMask, MovingImageMaskType);
/** Set/Get the fixed image mask. */
itkSetObjectMacro(FixedImageMask1, FixedImageMaskType);
itkSetObjectMacro(FixedImageMask2, FixedImageMaskType);
itkGetConstObjectMacro(FixedImageMask1, FixedImageMaskType);
itkGetConstObjectMacro(FixedImageMask2, FixedImageMaskType);
/** Set/Get gradient computation. */
itkSetMacro(ComputeGradient, bool);
itkGetConstReferenceMacro(ComputeGradient, bool);
itkBooleanMacro(ComputeGradient);
itkSetMacro(MaxTranslation, double);
itkGetConstReferenceMacro(MaxTranslation, double);
/** Get Gradient Image. */
itkGetConstObjectMacro(GradientImage, GradientImageType);
/** Set the parameters defining the Transform. */
bool
SetTransformParameters(const ParametersType& parameters) const;
/** Return the number of parameters required by the Transform */
unsigned int GetNumberOfParameters() const;
itk::Optimizer::ScalesType GetWeightings() const;
ParametersType GetParameters() const;
/** Initialize the Metric by making sure that all the components
* are present and plugged together correctly */
virtual void
Initialize();
protected:
gTwoImageToOneImageMetric();
~gTwoImageToOneImageMetric() override = default;
void
PrintSelf(std::ostream& os, Indent indent) const override;
constexpr static unsigned int Dimension = 3;
mutable unsigned long m_NumberOfPixelsCounted;
mutable double m_bestMeassure;
FixedImageConstPointer m_FixedImage1;
FixedImageConstPointer m_FixedImage2;
MovingImageConstPointer m_MovingImage;
mutable TransformPointer m_Transform1;
mutable TransformPointer m_Transform2;
mutable TransformPointer m_IsocTransf;
InterpolatorPointer m_Interpolator1;
InterpolatorPointer m_Interpolator2;
bool m_ComputeGradient;
GradientImagePointer m_GradientImage;
using OutputPixelType = unsigned char;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
mutable FixedImageMaskPointer m_FixedImageMask1;
mutable FixedImageMaskPointer m_FixedImageMask2;
mutable MovingImageMaskPointer m_MovingImageMask;
R23MetaInformation::Pointer
m_TransformMetaInfo;
InternalTransformMetaInformation::Pointer
m_internalMeta1,
m_internalMeta2;
ChangeInformationFilterPointer
m_Filter1,
m_Filter2;
double m_MaxTranslation;
/* Writes the image to the disk. If path is empty it will be saved to /img/out/*/
void WriteImage(MovingImageConstPointer img, std::string fileName, std::string path) const;
private:
FixedImageRegionType m_FixedImageRegion1;
FixedImageRegionType m_FixedImageRegion2;
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkgTwoImageToOneImageMetric.hxx"
#endif
#endif

476
reg23Topograms/itkDTRrecon/autoreg/itkgTwoImageToOneImageMetric.hxx Normal file
View File

@ -0,0 +1,476 @@
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkgTwoImageToOneImageMetric_hxx
#define itkgTwoImageToOneImageMetric_hxx
#include "itkgTwoImageToOneImageMetric.h"
#include <limits>
namespace itk {
template <typename TFixedImage, typename TMovingImage>
gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::gTwoImageToOneImageMetric()
{
m_FixedImage1 = nullptr; // has to be provided by the user.
m_FixedImage2 = nullptr; // has to be provided by the user.
m_MovingImage = nullptr; // has to be provided by the user.
m_Transform1 = nullptr; // has to be provided by the user.
m_Transform2 = nullptr; // has to be provided by the user.
m_TransformMetaInfo = nullptr; // has to be provided by the user.
m_internalMeta1 = nullptr;
m_internalMeta2 = nullptr;
m_Interpolator1 = nullptr; // has to be provided by the user.
m_Interpolator2 = nullptr; // has to be provided by the user.
m_GradientImage = nullptr; // will receive the output of the filter;
m_ComputeGradient = true; // metric computes gradient by default
m_NumberOfPixelsCounted = 0; // initialize to zero
m_bestMeassure = std::numeric_limits<short>::max();
m_MaxTranslation = std::numeric_limits<short>::max();
}
/*
* Set the parameters that define a unique transform
*/
template <typename TFixedImage, typename TMovingImage>
bool gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::SetTransformParameters(const ParametersType& parameters) const
{
if (!m_Transform1) {
itkExceptionMacro(<< "Transform 1 has not been assigned");
}
if (!m_Transform2) {
itkExceptionMacro(<< "Transform 2 has not been assigned");
}
if (!m_internalMeta1){
itkExceptionMacro(<< "m_internalMeta1 has not been assigned");
}
if (!m_internalMeta2){
itkExceptionMacro(<< "m_internalMeta2 has not been assigned");
}
if (!m_IsocTransf){
itkExceptionMacro(<< "m_IsocTransf has not been assigned");
}
//auto transformParameters = m_Transform1->GetParameters();
auto transformParameters = m_IsocTransf->GetParameters();
double TranslationAlongX = transformParameters[3];
double TranslationAlongY = transformParameters[4];
double TranslationAlongZ = transformParameters[5];
double RotationAlongX = transformParameters[0];
double RotationAlongY = transformParameters[1];
double RotationAlongZ = transformParameters[2];
//std::cout << "m_IsocTransf PARS: "<< transformParameters <<std::endl;
switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {
case X_ONLY:
TranslationAlongX = parameters[0];
break;
case Y_ONLY:
TranslationAlongY = parameters[0];
break;
case Z_ONLY:
TranslationAlongZ = parameters[0];
break;
case THREE_DOF:
TranslationAlongX = parameters[0];
TranslationAlongY = parameters[1];
TranslationAlongZ = parameters[2];
break;
case ROTATION_ONLY:
RotationAlongX = parameters[0];
RotationAlongY = parameters[1];
RotationAlongZ = parameters[2];
break;
case SIX_DOF:
TranslationAlongX = parameters[0];
TranslationAlongY = parameters[1];
TranslationAlongZ = parameters[2];
RotationAlongX = parameters[3];
RotationAlongY = parameters[4];
RotationAlongZ = parameters[5];
break;
default:
itkExceptionMacro(<< "Unsupported degree of freedeom");
break;
}
// std::cout << "New Transform Parameters ISOC IEC = " << std::endl;
// std::cout << " Translation X = " << TranslationAlongX << " mm" << std::endl;
// std::cout << " Translation Y = " << TranslationAlongY << " mm" << std::endl;
// std::cout << " Translation Z = " << TranslationAlongZ << " mm" << std::endl;
// std::cout << " Rotation Along X = " << RotationAlongX /*/ dtr*/ << " deg" << std::endl;
// std::cout << " Rotation Along Y = " << RotationAlongY /*/ dtr*/ << " deg" << std::endl;
// std::cout << " Rotation Along Z = " << RotationAlongZ /*/ dtr*/ << " deg" << std::endl;
// std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
// individual R and T components of the IsoC transform
ImageType3D::PointType
pT;
pT[0] = TranslationAlongX;
pT[1] = TranslationAlongY;
pT[2] = TranslationAlongZ;
ImageType3D::PointType
pR;
pR[0] = RotationAlongX;
pR[1] = RotationAlongY;
pR[2] = RotationAlongZ;
// check if we are within tolerance
bool transformValid =
(std::abs(pT[0]) < m_MaxTranslation) &&
(std::abs(pT[1]) < m_MaxTranslation) &&
(std::abs(pT[2]) < m_MaxTranslation);
if (transformValid) {
/*
Here we calculate the transform for eack projection corresponding to the isocentric
one we are optimizing.
This calculation takes into account calibrated center of projection for each view.
*/
// Calculate internal transform1
TransformType::Pointer CurrTransform;
CurrTransform = CalculateInternalTransformV3(
pT,
pR,
m_internalMeta1->GetpCalProjCenter(),
m_internalMeta1->GetpRTIsocenter(),
m_internalMeta1->GetIECtoLPSDirs());
m_Transform1->SetIdentity();
m_Transform1->SetParameters(
CurrTransform->GetParameters());
m_Transform1->SetCenter(
m_internalMeta1->GetpCalProjCenter());
// std::cout<<"----- itkgTwoImageToOne SetTransformParameters -----"<<std::endl;
// std::cout << "m_Transform1 :: Center: "<< m_Transform1->GetCenter()<<std::endl;
// std::cout << "m_Transform1 :: Pars: "<< m_Transform1->GetParameters()<<std::endl;
// std::cout<< "--------------------------------" <<std::endl;
// Calculate internal transform2
CurrTransform = CalculateInternalTransformV3(
pT,
pR,
m_internalMeta2->GetpCalProjCenter(),
m_internalMeta2->GetpRTIsocenter(),
m_internalMeta2->GetIECtoLPSDirs());
m_Transform2->SetIdentity();
m_Transform2->SetParameters(
CurrTransform->GetParameters());
m_Transform2->SetCenter(
m_internalMeta2->GetpCalProjCenter());
} else {
std::cout << " Transform invalid, out of range!" << std::endl;
std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
}
return transformValid;
}
template <typename TFixedImage, typename TMovingImage>
void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::Initialize()
{
if (!m_IsocTransf) {
itkExceptionMacro(<< "Transform is not present");
}
if (!m_Transform1) {
itkExceptionMacro(<< "Transform1 is not present");
}
if (!m_Transform2) {
itkExceptionMacro(<< "Transform2 is not present");
}
if(!m_internalMeta1) {
itkExceptionMacro(<< "m_internalMeta1 is not present");
}
if(!m_internalMeta2) {
itkExceptionMacro(<< "m_internalMeta2 is not present");
}
if (!m_Interpolator1) {
itkExceptionMacro(<< "Interpolator1 is not present");
}
if (!m_Interpolator2) {
itkExceptionMacro(<< "Interpolator2 is not present");
}
if (!m_FixedImage1) {
itkExceptionMacro(<< "FixedImage1 is not present");
}
if (!m_FixedImage2) {
itkExceptionMacro(<< "FixedImage2 is not present");
}
if (m_FixedImageRegion1.GetNumberOfPixels() == 0) {
itkExceptionMacro(<< "FixedImageRegion1 is empty");
}
if (m_FixedImageRegion2.GetNumberOfPixels() == 0) {
itkExceptionMacro(<< "FixedImageRegion2 is empty");
}
// If the image is provided by a source, update the source.
if (m_MovingImage->GetSource()) {
m_MovingImage->GetSource()->Update();
}
// If the image is provided by a source, update the source.
if (m_FixedImage1->GetSource()) {
m_FixedImage1->GetSource()->Update();
}
if (m_FixedImage2->GetSource()) {
m_FixedImage2->GetSource()->Update();
}
//std::cout<<"FixedImageRegion1 " <<this->GetFixedImageRegion1() <<std::endl;
//std::cout<<"m_FixedImage1 buffered " <<this->m_FixedImage1->GetBufferedRegion() <<std::endl;
//std::cout<<"Moving buffered " <<this->m_Filter1->GetOutput()->GetBufferedRegion() <<std::endl;
//std::cout<<"GetFixedImage1 buffered " <<this->GetFixedImage1()->GetBufferedRegion() <<std::endl;
//std::cout<<"Filter1 buffered " <<this->m_Filter1->GetOutput()->GetBufferedRegion() <<std::endl;
// Make sure the FixedImageRegion is within the FixedImage buffered region
if (!m_FixedImageRegion1.Crop(m_FixedImage1->GetBufferedRegion())) {
itkExceptionMacro(<< "FixedImageRegion1 does not overlap the fixed image buffered region");
}
if (!m_FixedImageRegion2.Crop(m_FixedImage2->GetBufferedRegion())) {
itkExceptionMacro(<< "FixedImageRegion2 does not overlap the fixed image buffered region");
}
this->SetTransformParameters(this->GetParameters());
if (m_ComputeGradient) {
GradientImageFilterPointer gradientFilter = GradientImageFilterType::New();
gradientFilter->SetInput(m_MovingImage);
const typename MovingImageType::SpacingType& spacing = m_MovingImage->GetSpacing();
double maximumSpacing = 0.0;
for (unsigned int i = 0; i < MovingImageDimension; i++) {
if (spacing[i] > maximumSpacing) {
maximumSpacing = spacing[i];
}
}
gradientFilter->SetSigma(maximumSpacing);
gradientFilter->SetNormalizeAcrossScale(true);
gradientFilter->Update();
m_GradientImage = gradientFilter->GetOutput();
}
// If there are any observers on the metric, call them to give the
// user code a chance to set parameters on the metric
this->InvokeEvent(InitializeEvent());
}
template <typename TFixedImage, typename TMovingImage>
unsigned int gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetNumberOfParameters() const
{
if (!m_TransformMetaInfo) {
itkExceptionMacro(<< "TransformMetaInfo is missing");
}
switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {
case X_ONLY:
case Y_ONLY:
case Z_ONLY:
return 1;
case THREE_DOF:
case ROTATION_ONLY:
return 3;
case SIX_DOF:
return 6;
default:
itkExceptionMacro(<< "Unsupported degree of freedeom");
break;
}
return m_IsocTransf->GetNumberOfParameters();
}
template <typename TFixedImage, typename TMovingImage>
itk::Optimizer::ScalesType gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetWeightings() const
{
itk::Optimizer::ScalesType weightings(this->GetNumberOfParameters());
switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {
case X_ONLY:
case Y_ONLY:
case Z_ONLY:
weightings[0] = 1.0;
break;
;
case THREE_DOF:
weightings[0] = 1.;
weightings[1] = 1.;
weightings[2] = 1.;
break;
case ROTATION_ONLY:
weightings[0] = 1. / dtr;
weightings[1] = 1. / dtr;
weightings[2] = 1. / dtr;
break;
case SIX_DOF:
//G: ITHINKTHOSEARE FLIPPED!!
weightings[0] = 1.;
weightings[1] = 1.;
weightings[2] = 1.;
weightings[3] = 1. / dtr;
weightings[4] = 1. / dtr;
weightings[5] = 1. / dtr;
break;
default:
itkExceptionMacro(<< "Unsupported degree of freedeom");
break;
}
return weightings;
}
template <typename TFixedImage, typename TMovingImage>
typename gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::
ParametersType gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetParameters() const
{
// ParametersType parametersTransform = m_Transform1->GetParameters(); //angleX, angleY, angleZ, transX, transY, transZ
ParametersType parametersTransform = m_IsocTransf->GetParameters(); //angleX, angleY, angleZ, transX, transY, transZ
ParametersType parameters(this->GetNumberOfParameters());
switch (m_TransformMetaInfo->GetDegreeOfFreedom()) {
case X_ONLY:
parameters[0] = parametersTransform[3];
break;
case Y_ONLY:
parameters[0] = parametersTransform[4];
break;
case Z_ONLY:
parameters[0] = parametersTransform[5];
break;
case THREE_DOF:
parameters[0] = parametersTransform[3];
parameters[1] = parametersTransform[4];
parameters[2] = parametersTransform[5];
break;
case ROTATION_ONLY:
parameters[0] = parametersTransform[0];
parameters[1] = parametersTransform[1];
parameters[2] = parametersTransform[2];
break;
case SIX_DOF:
parameters[0] = parametersTransform[3];
parameters[1] = parametersTransform[4];
parameters[2] = parametersTransform[5];
parameters[3] = parametersTransform[0];
parameters[4] = parametersTransform[1];
parameters[5] = parametersTransform[2];
break;
default:
itkExceptionMacro(<< "Unsupported degree of freedeom");
break;
}
return parameters;
}
template <typename TFixedImage, typename TMovingImage>
void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::PrintSelf(std::ostream& os, Indent indent) const
{
Superclass::PrintSelf(os, indent);
os << indent << "ComputeGradient: " << static_cast<typename NumericTraits<bool>::PrintType>(m_ComputeGradient)
<< std::endl;
os << indent << "Moving Image: " << m_MovingImage.GetPointer() << std::endl;
os << indent << "Fixed Image 1: " << m_FixedImage1.GetPointer() << std::endl;
os << indent << "Fixed Image 2: " << m_FixedImage2.GetPointer() << std::endl;
os << indent << "Gradient Image: " << m_GradientImage.GetPointer() << std::endl;
os << indent << "Transform1: " << m_Transform1.GetPointer() << std::endl;
os << indent << "Transform2: " << m_Transform2.GetPointer() << std::endl;
os << indent << "Interpolator 1: " << m_Interpolator1.GetPointer() << std::endl;
os << indent << "Interpolator 2: " << m_Interpolator2.GetPointer() << std::endl;
os << indent << "Filter 1: " << m_Filter1.GetPointer() << std::endl;
os << indent << "Filter 2: " << m_Filter2.GetPointer() << std::endl;
os << indent << "Transform MetaInfoe: " << m_TransformMetaInfo.GetPointer() << std::endl;
os << indent << "FixedImageRegion 1: " << m_FixedImageRegion1 << std::endl;
os << indent << "FixedImageRegion 2: " << m_FixedImageRegion2 << std::endl;
os << indent << "Moving Image Mask: " << m_MovingImageMask.GetPointer() << std::endl;
os << indent << "Fixed Image Mask 1: " << m_FixedImageMask1.GetPointer() << std::endl;
os << indent << "Fixed Image Mask 2: " << m_FixedImageMask2.GetPointer() << std::endl;
os << indent << "Number of Pixels Counted: " << m_NumberOfPixelsCounted << std::endl;
}
template <typename TFixedImage, typename TMovingImage>
void gTwoImageToOneImageMetric<TFixedImage, TMovingImage>::WriteImage(MovingImageConstPointer img, std::string fileName, std::string path) const
{
using RescaleFilterType = itk::RescaleIntensityImageFilter<MovingImageType, OutputImageType>;
using WriterType = itk::ImageFileWriter<OutputImageType>;
if (img) {
// Rescale the intensity of the projection images to 0-255 for output.
typename RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
rescaler1->SetOutputMinimum(0);
rescaler1->SetOutputMaximum(255);
rescaler1->SetInput(img);
WriterType::Pointer writer1 = WriterType::New();
if (fileName.empty()) {
char result[PATH_MAX];
ssize_t count = readlink("/proc/self/exe", result, PATH_MAX);
if (count != -1) {
path = dirname(result);
path += "/img/out/";
}
}
if (fileName.empty()) {
fileName = path;
fileName += "2D1.tif";
}
writer1->SetFileName(fileName);
writer1->SetInput(rescaler1->GetOutput());
try {
std::cout << "Writing image: " << fileName << std::endl;
writer1->Update();
} catch (itk::ExceptionObject& err) {
std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
}
}
} // end namespace itk
#endif

1505
reg23Topograms/itkDTRrecon/itkImageProcessor.cpp
View File

File diff suppressed because it is too large Load Diff

150
reg23Topograms/itkDTRrecon/itkImageProcessor.h
View File

@ -38,20 +38,18 @@ gfattori 08.11.2021
#include "gdcmGlobal.h"
#include "vtkImageData.h"
#include "vtkTransform.h"
#include "vtkSmartPointer.h"
#include "itkQtIterationUpdate.h"
#include "itkImageProcessorHelpers.h"
#include "itkReg23.h"
namespace itk
{
/* reference string required for comparison with tag values */
static const char *ImageOrientationStrings[] = {
"NotDefined",
"HFS",
"FFS"
};
class ITK_EXPORT itkImageProcessor : public itk::Object
{
@ -71,13 +69,50 @@ public:
/** Run-time type information (and related methods). */
itkTypeMacro(itkImageProcessor, Object);
CommandIterationUpdate::Pointer
optimizerObserver;
ExhaustiveCommandIterationUpdate::Pointer
ExhaustiveOptimizerObserver;
/** Sets the degree of freedom for optimzation*/
tDegreeOfFreedomEnum GetDegreeOfFreedom();
/** Returns user components only of the autoReg
* */
Optimizer::ParametersType
GetFinalR23Parameters();
/** Define the region to be used as ROI on fixed images */
void SetRegionFixed1(int,int,int,int);
void SetRegionFixed2(int,int,int,int);
void ResetROIRegions();
void InizializeRegistration(double stepLength,
double maxTranslation,
eDegreeOfFreedomType dof);
/** Maximum number of iterations for the optimizer */
void SetMaxNumberOfIterations(int);
itkGetMacro(R23, itkReg23::Pointer);
/** Set Optimizer type */
void SetOptimizer(std::string);
/** Set Metric type */
void SetMetric(std::string);
/** Set number of logic CPU to be made available to interpolators*/
void SetNumberOfWorkingUnits(int iN);
/** 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);
@ -96,8 +131,14 @@ public:
double GetSCD1();
double GetSCD2();
void SetCustom_ProjectionCenterOffsetFixedAxes_IEC(double,double,double,double,double,double);
/** 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();
void SetCustom_ProjectionCenterOffsetFixedAxes_IEC(double,double,double,double,double,double);
void SetCustom_ProjectionCenterFixedAxes_IEC(double,double,double);
/** Intensity threshold used for image projection,
@ -107,14 +148,49 @@ public:
/** Custom settings of the projection images */
void SetCustom_2Dres(double,double,double,double);
void SetCustom_2Dsize(int,int,int,int);
/** Custom transform to map IEC of imaging device into IEC-Support
* Tx Ty Tz [mm] Rx Ry Rz [deg]
*/
void SetCustom_ImportTransform(double, double, double,
double, double, double);
/** Should rotations be used when computing import offset and hidden transform? */
void SetUseRotationsForImportOffset(bool bVal);
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 */
/** Get the complete transform, likely for export to SRO.
* This includes user and hidden contributions.
* Transform center has to be specified, Zero would result
* in Isocentric transform with center at the RT Iso.
*/
TransformType::Pointer
GetCompleteIsocentricTransform(/*ImageType3D::PointType TransformCenter*/);
/** Get only the User component of the transform
* as parameters' list. Import offset and hidden not included.
* This can be use to update UI transform display
* ParametersType :: Rx Ry Rz Tx Ty Tz
*/
Optimizer::ParametersType
GetUserIsocentricTransform();
/** Return useful meta info for external use
* Consumer should check for nullptr...
*/
const CTVolumeImageMetaInformation::Pointer
GetCTMetaInfo();
const RTGeometryMetaInformation::Pointer
GetRTMetaInfo();
const DRTProjectionGeometryImageMetaInformation::Pointer
GetDRTGeoMetaInfo();
/** Initialize projection geometry */
void InitializeProjector();
@ -125,8 +201,6 @@ public:
/** Get Projection origin in LPS coordinates*/
const std::vector <double> GetRTImportOffset();
/** Get the Localizer intensity window and
* center for rendering */
double GetLocalizerDisplayWindowLevel(int );
@ -150,13 +224,12 @@ public:
void Write2DImages();
protected:
/** Various pixel types */
using InternalPixelType = float;
using PixelType2D = short;
using PixelType3D = short;
using OutputPixelType = unsigned char;
using OutputPixelType = unsigned short;
using ImageType3D = itk::Image<PixelType3D, Dimension>;
using InternalImageType = itk::Image<InternalPixelType, Dimension>;
@ -181,11 +254,12 @@ private:
/** 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>;
using RoiForRegistration = itk::ImageRegion<Dimension>;
/** Image reader types */
using ImageReaderType2D = itk::ImageFileReader<ImageType2D>;
using ImageReaderType3D = itk::ImageFileReader<ImageType3D>;
@ -202,23 +276,10 @@ private:
/** ITK to VTK filter */
using ITKtoVTKFilterType = itk::ImageToVTKImageFilter<OutputImageType>;
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
);
TransformType::Pointer
transform1,
transform2;
transform2,
IsocTransf;
InterpolatorType::Pointer
interpolator1,
@ -235,6 +296,8 @@ private:
imageDRT1In,
imageDRT2In;
DuplicatorType::Pointer
m_LATSourceDupli,
m_PASourceDupli,
@ -265,13 +328,6 @@ private:
);
TransformType::Pointer
MapTransformToNewOrigin(
ImageType3D::PointType m_COR,
ImageType3D::PointType m_Translations,
ImageType3D::PointType m_Rotations
);
double
CalcProjectionAngleLPS(
tPatOrientation pOrient,
@ -305,8 +361,6 @@ private:
* m_Projection1VTKTransform,
* m_Projection2VTKTransform;
/**
* Many meta containers
*/
@ -330,6 +384,26 @@ private:
TransformMetaInformation::Pointer
m_TransformMetaInfo;
R23MetaInformation::Pointer
m_r23MetaInfo;
InternalTransformMetaInformation::Pointer
m_InternalTransf1,
m_InternalTransf2;
RoiForRegistration
roiAutoReg1,
roiAutoReg2;
int
iNWUnits;
itk::itkReg23::Pointer
m_R23;
};

165
reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.cpp Normal file
View File

@ -0,0 +1,165 @@
#include "itkImageProcessorHelpers.h"
#include <gdcmReader.h>
#include <stdio.h>
namespace itk {
TransformType::Pointer
MapTransformToNewOrigin(
ImageType3D::PointType m_COR, // Center of rotation for the proj geometry. this is my new origin.
ImageType3D::PointType m_Translations,
ImageType3D::PointType m_Rotations
){
TransformType::Pointer InputTransform = TransformType::New();
InputTransform->SetComputeZYX(true);
InputTransform->SetIdentity();
TransformType::OutputVectorType translation;
translation[0] = m_Translations[0];
translation[1] = m_Translations[1];
translation[2] = m_Translations[2];
InputTransform->SetTranslation(translation);
const double dtr = (atan(1.0) * 4.0) / 180.0;
InputTransform->SetRotation(
dtr * m_Rotations[0],
dtr * m_Rotations[1],
dtr * m_Rotations[2]);
ImageType3D::PointType m_TransformOrigin;
m_TransformOrigin.Fill(0.);
InputTransform->SetCenter(
m_TransformOrigin );
ImageType3D::PointType NewOriginTranslations =
InputTransform->TransformPoint(m_COR);
ImageType3D::PointType DeltaNewOrigin =
NewOriginTranslations - m_COR;
TransformType::Pointer m_OutputTransform =
TransformType::New();
m_OutputTransform ->SetComputeZYX(true);
m_OutputTransform ->SetIdentity();
translation[0] = DeltaNewOrigin[0];
translation[1] = DeltaNewOrigin[1];
translation[2] = DeltaNewOrigin[2];
m_OutputTransform->SetTranslation(translation);
m_OutputTransform->SetRotation(
dtr * m_Rotations[0],
dtr * m_Rotations[1],
dtr * m_Rotations[2]);
m_OutputTransform->SetCenter(m_COR);
InputTransform = NULL;
return m_OutputTransform;
}
TransformType::Pointer
CalculateInternalTransformV3(
ImageType3D::PointType m_Translation, //IEC
ImageType3D::PointType m_Rotation, //IEC
ImageType3D::PointType m_CalibratedProjectionCenter, //LPS
ImageType3D::PointType m_RTIsocenter, //LPS
InternalImageType::DirectionType m_IECtoLPSDirections
)
{
//Convert all inputs into LPS
ImageType3D::PointType m_TLPS =
m_IECtoLPSDirections * m_Translation;
ImageType3D::PointType m_RLPS =
m_IECtoLPSDirections * m_Rotation;
// Map offset to the projection center
TransformType::Pointer m_outputTransform =
MapTransformToNewOrigin (
m_CalibratedProjectionCenter - m_RTIsocenter,
m_TLPS,
m_RLPS
);
m_outputTransform->SetCenter(m_CalibratedProjectionCenter);
return m_outputTransform;
}
const char* queryStationName(const char* pcFName){
static std::string buffer;
/* Check if we can open the file */
gdcm::Reader R;
R.SetFileName(pcFName);
if(!R.Read())
{std::cerr<<"ERROR: cannot read file: "<<pcFName<<std::endl;
return "\n";
}
const gdcm::File &file = R.GetFile();
const gdcm::DataSet &ds = file.GetDataSet();
std::string s (gGetStringValueFromTag(gdcm::Tag(0x0008,0x1010), ds));
std::copy(s.rbegin(), s.rend(), std::back_inserter(buffer));
std::reverse(buffer.begin(),buffer.end());
return buffer.c_str();
}
int query2DimageReconstructionDiameter(const char *pcFName){
/* Check if we can open the file */
gdcm::Reader R;
R.SetFileName(pcFName);
if(!R.Read())
{std::cerr<<"ERROR: cannot read file: "<<pcFName<<std::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,0x1100), ds));
return std::atoi(sTmpString);
}
//FUNCTION DECLARATION NECESSARY FOR COPY/PASTE FROM vtkGDCMImageReader
// const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds);
const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds)
{
static std::string buffer;
buffer = ""; // cleanup previous call
if( ds.FindDataElement( t ) )
{
const gdcm::DataElement& de = ds.GetDataElement( t );
const gdcm::ByteValue *bv = de.GetByteValue();
if( bv ) // Can be Type 2
{
buffer = std::string( bv->GetPointer(), bv->GetLength() );
// Will be padded with at least one \0
}
}
// Since return is a const char* the very first \0 will be considered
return buffer.c_str();
}
}

55
reg23Topograms/itkDTRrecon/itkImageProcessorHelpers.h Normal file
View File

@ -0,0 +1,55 @@
#ifndef ITKIMAGEPROCESSORHELPERS_H
#define ITKIMAGEPROCESSORHELPERS_H
#include "itkEuler3DTransform.h"
#include "itkImage.h"
#include "itkEuler3DTransform.h"
#include "itkCommand.h"
#include "itkObject.h"
#include "itkObjectFactory.h"
#include <gdcmTag.h>
#include <gdcmDataSet.h>
namespace itk
{
using TransformType = itk::Euler3DTransform<double>;
constexpr static unsigned int Dimension = 3;
using PixelType3D = short;
using InternalPixelType = float;
using ImageType3D = itk::Image<PixelType3D, Dimension>;
using InternalImageType = itk::Image<InternalPixelType, Dimension>;
TransformType::Pointer
MapTransformToNewOrigin(
ImageType3D::PointType m_COR,
ImageType3D::PointType m_Translations,
ImageType3D::PointType m_Rotations
);
TransformType::Pointer
CalculateInternalTransformV3(
ImageType3D::PointType m_Translation,
ImageType3D::PointType m_Rotation,
ImageType3D::PointType m_CalibratedProjectionCenter,
ImageType3D::PointType m_RTIsocenter,
InternalImageType::DirectionType m_IECtoLPSDirections
);
const char *gGetStringValueFromTag(const gdcm::Tag& t, const gdcm::DataSet& ds);
int query2DimageReconstructionDiameter(const char*);
const char* queryStationName(const char*);
}
#endif

192
reg23Topograms/itkDTRrecon/itkQtIterationUpdate.h Normal file
View File

@ -0,0 +1,192 @@
#ifndef ITKQTITERATIONUPDATE_H
#define ITKQTITERATIONUPDATE_H
#include "itkCommand.h"
#include <QObject>
#include "itkPowellOptimizer.h"
#include "itkAmoebaOptimizer.h"
#include "itkExhaustiveOptimizer.h"
#include "itkTwoProjectionImageRegistrationMethod.h"
class QObjectIterationUpdate : public QObject{
Q_OBJECT
public:
QObjectIterationUpdate(){
bAbortProcessCommand = false;
};
bool getAbortFlag(){
return bAbortProcessCommand;
};
void setAbortFlag(bool bVal){
bAbortProcessCommand = bVal;
};
public slots:
void onAbortProcessCommand(){
bAbortProcessCommand = true;
};
void onIteration(double dI,double dX,double dY,double dZ){
emit
sendRegistrationProgress(dI,dX,dY,dZ);
}
private:
bool
bAbortProcessCommand;
signals:
void sendRegistrationProgress(double,double,double,double);
};
namespace itk
{
class CommandIterationUpdate : public itk::Command {
// TODO: Move to own files.
constexpr static unsigned int Dimension = 3;
public:
QObjectIterationUpdate *objIterUpdate;
using Self = CommandIterationUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(CommandIterationUpdate);
using InternalPixelType = float;
using InternalImageType = itk::Image<InternalPixelType, Dimension>;
using ProcesssType = typename itk::TwoProjectionImageRegistrationMethod<InternalImageType, InternalImageType>;
using ProcesssPointer = typename ProcesssType::Pointer;
/** Set/Get the Process. */
itkSetObjectMacro(Process, ProcesssType);
itkGetConstObjectMacro(Process, ProcesssType);
private:
protected:
CommandIterationUpdate() {
objIterUpdate = new QObjectIterationUpdate;
}
ProcesssPointer m_Process;
public:
using OptimizerType = itk::PowellOptimizer;
using OptimizerPointer = const OptimizerType*;
using AmoebaOptimizerType = itk::AmoebaOptimizer;
using AmoebaOptimizerPointer = const OptimizerType*;
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
{
if(objIterUpdate->getAbortFlag()){
std::cout << "AbortRequest" << std::endl;
objIterUpdate->setAbortFlag(false);
throw itk::ProcessAborted();
}
// std::cout << "Progress: " << this->m_Process->GetAbortGenerateData() << std::endl;
auto optimizer = dynamic_cast<OptimizerPointer>(object);
if (typeid(event) == typeid(itk::IterationEvent)) {
//Feedback from the optimizer executed at the end of every itteration
// currently just print the result into the cout. We might add
// functionality to register and emit signals to update the UI.
std::cout << "Iteration: " << optimizer->GetCurrentIteration() << std::endl;
auto oldprecision = std::cout.precision();
std::cout.precision(std::numeric_limits<double>::digits10 + 2);
std::cout << "Similarity: " << optimizer->GetValue() << std::endl;
std::cout.precision(oldprecision);
std::cout << "Position: " << optimizer->GetCurrentPosition() << std::endl;
// objIterUpdate->onIteration(
// optimizer->GetCurrentIteration()+1,
// optimizer->GetCurrentPosition()[0],
// optimizer->GetCurrentPosition()[2],
// -optimizer->GetCurrentPosition()[1]
// );
objIterUpdate->onIteration(
optimizer->GetCurrentIteration()+1,
optimizer->GetCurrentPosition()[0],
optimizer->GetCurrentPosition()[1],
optimizer->GetCurrentPosition()[2]
);
}
return;
}
};
class ExhaustiveCommandIterationUpdate : public itk::Command {
// TODO: Move to own files.
public:
using Self = ExhaustiveCommandIterationUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
std::ostream* os;
protected:
ExhaustiveCommandIterationUpdate() = default;
public:
using OptimizerType = itk::ExhaustiveOptimizer;
;
using OptimizerPointer = const OptimizerType*;
void set_stream(std::ostream& stream)
{
os = &stream;
}
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)) {
//crude LPS to IEC transform for HFS.
auto position = optimizer->GetCurrentPosition();
auto oldprecision = os->precision();
os->precision(std::numeric_limits<double>::digits10 + 2);
*os << optimizer->GetCurrentValue();
os->precision(oldprecision);
//*os << "\t" << position[0] << "\t" << position[2] << "\t" << -position[1] << std::endl;
*os << "\t" << position[0] << "\t" << position[1] << "\t" << position[2] << std::endl;
}
return;
}
};
}
#endif

393
reg23Topograms/itkDTRrecon/itkReg23.cpp Normal file
View File

@ -0,0 +1,393 @@
#include "itkReg23.h"
#include "itkTimeProbesCollectorBase.h"
namespace itk
{
itkReg23::itkReg23()
{
// Metrics
NCCmetric = MetricType::New();
MImetric = MIMetricType::New();
// Optimizers
PowellOptimizer = PowellOptimizerType::New();
AmoebaOptimizer = AmoebaOptimizerType::New();
ExhaustiveOptimizer = ExhaustiveOptimizerType::New();
// Registration
registration = RegistrationType::New();
//internalTransforms
Transform1 = TransformType::New();
Transform2 = TransformType::New();
IsocTransform = TransformType::New();
registration->SetTransform1(Transform1);
registration->SetTransform2(Transform2);
// to be provided by the user
m_r23Meta = nullptr;
m_Volume = nullptr;
m_PA = nullptr;
m_LAT = nullptr;
m_InternalTransf1 = nullptr;
m_InternalTransf2 = nullptr;
m_filter1 = nullptr;
m_filter2 = nullptr;
m_TransformMetaInfo = nullptr;
m_OptimizerObserver = nullptr;
// m_UseFullROI = true; // if the full image ROI shall be used
m_UseDumptoFile = false; //If each (!) optimzer result shall be dumpped into a file.
}
itkReg23::~itkReg23()
{
}
void itkReg23::PrintSelf(std::ostream& os, Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
void itkReg23::InitializeRegistration(
double stepLength,
double maxTranslation,
eDegreeOfFreedomType dof)
{
std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
if (!m_PA) {
itkExceptionMacro(<< "PA topogram not present");
}
if (!m_LAT) {
itkExceptionMacro(<< "LAT topogram not present");
}
if (!m_Volume) {
itkExceptionMacro(<< "CT data not present");
}
if (!m_r23Meta) {
itkExceptionMacro(<< "r23Meta data not present");
}
if (!m_OptimizerObserver) {
itkExceptionMacro(<< "m_OptimizerObserver data not present");
}
if (!m_InternalTransf1) {
itkExceptionMacro(<< "m_InternalTransf1 data not present");
}
if (!m_InternalTransf2) {
itkExceptionMacro(<< "m_InternalTransf2 data not present");
}
if (!m_filter1) {
itkExceptionMacro(<< "m_filter1 data not present");
}
if (!m_filter2) {
itkExceptionMacro(<< "m_filter2 data not present");
}
if (!m_TransformMetaInfo) {
itkExceptionMacro(<< "m_TransformMetaInfo data not present");
}
if(!m_interpolator1){
itkExceptionMacro(<< "m_interpolator1 data not present");
}
if(!m_interpolator2){
itkExceptionMacro(<< "m_interpolator2 data not present");
}
AmoebaOptimizerType::ParametersType simplexDelta(3);
ExhaustiveOptimizerType::StepsType steps(3);
const int numberOfSteps = 25; //In each direction. Total number of steps is ((2*numberOfSteps+1))^3. For 25 -> 132651.
//const double stepLength = 0.1;
switch (m_r23Meta->GetOptimizerType()) {
case tOptimizerTypeEnum::POWELL:
std::cout<< "Using POWELL Optimizer" <<std::endl;
PowellOptimizer->SetMaximize(false); // for NCC and MI
PowellOptimizer->SetMaximumIteration(m_r23Meta->GetMaxIterations());
PowellOptimizer->SetMaximumLineIteration(4); // for Powell's method
PowellOptimizer->SetStepLength(stepLength);
PowellOptimizer->SetStepTolerance(0.01);
PowellOptimizer->SetValueTolerance(0.000002);
PowellOptimizer->RemoveAllObservers();
PowellOptimizer->AddObserver(itk::IterationEvent(), m_OptimizerObserver);
registration->SetOptimizer(PowellOptimizer);
break;
case tOptimizerTypeEnum::AMOEBA:
std::cout<< "Using AMOEBA Optimizer" <<std::endl;
AmoebaOptimizer->SetMinimize(true);
AmoebaOptimizer->SetMaximumNumberOfIterations(m_r23Meta->GetMaxIterations());
AmoebaOptimizer->SetParametersConvergenceTolerance(0.1);
AmoebaOptimizer->SetFunctionConvergenceTolerance(0.000002);
AmoebaOptimizer->SetAutomaticInitialSimplex(false);
//Initial size of the simplex (otherwise it is a very small number and optimizer stops immeditaly)
// 2 mm / 2 degrees seems to be a good value which performs well.
if (GetNumberOfDOF(dof) == 0) {
itkExceptionMacro(<< "Unkown or unsupported degree of freedom");
}
for (int i = 0; i < GetNumberOfDOF(dof); i++) {
simplexDelta[i] = 2.0;
}
AmoebaOptimizer->SetInitialSimplexDelta(simplexDelta);
AmoebaOptimizer->RemoveAllObservers();
AmoebaOptimizer->AddObserver(itk::IterationEvent(), m_OptimizerObserver);
registration->SetOptimizer(AmoebaOptimizer);
break;
case tOptimizerTypeEnum::EXHAUSTIVE:
std::cout<< "Using Extensive Optimizer" <<std::endl;
steps[0] = numberOfSteps;
steps[1] = numberOfSteps;
steps[2] = numberOfSteps;
ExhaustiveOptimizer->SetNumberOfSteps(steps);
ExhaustiveOptimizer->SetStepLength(stepLength);
registration->SetOptimizer(ExhaustiveOptimizer);
break;
default:
break;
}
switch (m_r23Meta->GetMetricType()) {
case tMetricTypeEnum::MI:
std::cout<< "Using MI Metric" <<std::endl;
registration->SetMetric(MImetric);
MImetric->ComputeGradientOff();
MImetric->SetMaxTranslation(maxTranslation);
break;
case tMetricTypeEnum::NCC:
std::cout<< "Using NCC Metric" <<std::endl;
registration->SetMetric(NCCmetric);
NCCmetric->ComputeGradientOff();
NCCmetric->SetSubtractMean(true);
NCCmetric->SetMaxTranslation(maxTranslation);
break;
default:
break;
}
registration->SetFixedImage1(m_PA);
registration->SetFixedImage2(m_LAT);
registration->SetMovingImage(m_Volume);
registration->SetFixedImageRegion1(m_roiAutoReg1);
registration->SetFixedImageRegion2(m_roiAutoReg2);
registration->SetTransformMetaInfo(m_r23Meta);
registration->SetinternalMeta1(m_InternalTransf1);
registration->SetinternalMeta2(m_InternalTransf2);
registration->SetInterpolator1(m_interpolator1);
registration->SetInterpolator2(m_interpolator2);
registration->SetFilter1(m_filter1);
registration->SetFilter2(m_filter2);
IsocTransform = TransformType::New();
IsocTransform->SetComputeZYX(true);
IsocTransform->SetIdentity();
IsocTransform->SetRotation(
m_TransformMetaInfo->GetR()[0],
m_TransformMetaInfo->GetR()[1],
m_TransformMetaInfo->GetR()[2]
);
TransformType::OutputVectorType TranslV;
TranslV[0] = m_TransformMetaInfo->GetT()[0];
TranslV[1] = m_TransformMetaInfo->GetT()[1];
TranslV[2] = m_TransformMetaInfo->GetT()[2];
IsocTransform->SetTranslation(TranslV);
registration->SetIsocIECTransform(IsocTransform);
// if (verbose) {
// registration->DebugOn();
// registration->Print(std::cout);
// }
std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
}
int itkReg23::StartRegistration(std::string extraInfo)
{
std::cout << "*" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
// TODO: Check if the registartion pipeline has been initialized
using ParametersType = RegistrationType::ParametersType;
//auto startParameters = transform1->GetParameters();
time_t t = time(0); // get time now
struct tm* now = localtime(&t);
bool useDumpToFile = false;
char buffer[255];
strftime(buffer, 255, "test-%F-%H%M%S.txt", now);
std::fstream fs;
if (useDumpToFile) {
fs.open(buffer, std::fstream::out);
fs << extraInfo;
fs << "Value\tX\tY\tZ " << std::endl;
// if (r23Meta->GetOptimizerType() == tOptimizerTypeEnum::EXHAUSTIVE) {
// ExhaustiveOptimizerObserver->set_stream(fs);
// }
}
// Start the registration
// ~~~~~~~~~~~~~~~~~~~~~~
// Create a timer to record calculation time.
itk::TimeProbesCollectorBase timer;
if (verbose) {
std::cout << "Starting the registration now..." << std::endl;
}
try {
// timer.Start("Registration");
// Start the registration.
registration->StartRegistration();
// timer.Stop("Registration");
} catch (itk::ExceptionObject& err) {
registration->ResetPipeline();
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return -1;
}
fs.close();
itk::Optimizer::ParametersType finalParameters =
this->GetCurrentPosition();
std::cout<<" FinalPars: "<< finalParameters <<std::endl;
std::cout << "#" << __COMPACT_PRETTY_FUNCTION__ << std::endl;
return 0;
}
double itkReg23::GetOptimizerValue()
{
return m_OptmizerValue;
}
/** Get the cost function value for the current transform*/
Optimizer::ParametersType itkReg23::GetCurrentPosition(){
itk::Optimizer::ParametersType finalParameters;
switch (m_r23Meta->GetOptimizerType()) {
case tOptimizerTypeEnum::POWELL:
finalParameters = PowellOptimizer->GetCurrentPosition();
m_OptmizerValue = PowellOptimizer->GetValue();
break;
case tOptimizerTypeEnum::AMOEBA:
finalParameters = AmoebaOptimizer->GetCurrentPosition();
m_OptmizerValue = AmoebaOptimizer->GetValue();
break;
case tOptimizerTypeEnum::EXHAUSTIVE:
finalParameters = ExhaustiveOptimizer->GetCurrentPosition();
break;
default:
break;
}
return finalParameters;
}
double itkReg23::GetCurrentMetricValue()
{
switch (m_r23Meta->GetMetricType()) {
case tMetricTypeEnum::MI:
if(!MImetric){
itkExceptionMacro(<< "MImetric not present");
return -1.;
} else {
return MImetric->GetValue();
}
break;
case tMetricTypeEnum::NCC:
if(!NCCmetric){
itkExceptionMacro(<< "NCCmetric not present");
return -1.;
} else {
return NCCmetric->GetValue();
}
break;
default:
break;
}
return -1.;
}
}

197
reg23Topograms/itkDTRrecon/itkReg23.h Normal file
View File

@ -0,0 +1,197 @@
#ifndef ITKREG23_H
#define ITKREG23_h
#include "itkCommand.h"
#include "itkExhaustiveOptimizer.h"
#include "itkMutualInformationTwoImageToOneImageMetric.h"
#include "itkNormalizedCorrelationTwoImageToOneImageMetric.h"
#include "itkPowellOptimizer.h"
#include "itkTwoProjectionImageRegistrationMethod.h"
#include "itkAmoebaOptimizer.h"
#include "itkCommand.h"
#include "itkObject.h"
#include "itkObjectFactory.h"
#include "itkSmartPointer.h"
#include "itkImageProcessorHelpers.h"
#include "itkQtIterationUpdate.h"
namespace itk{
class ITK_EXPORT itkReg23 : public itk::Object
{
public:
/** Standard "Self" typedef. */
typedef itkReg23 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(itkReg23, Object);
using ChangeInformationFilterType = itk::ChangeInformationImageFilter<InternalImageType>;
using RoiForRegistration = itk::ImageRegion<Dimension>;
using InterpolatorType = itk::gSiddonJacobsRayCastInterpolateImageFunction<InternalImageType, double>;
/* ---- User provided */
itkSetMacro(OptimizerObserver, CommandIterationUpdate::Pointer);
itkGetMacro(OptimizerObserver, CommandIterationUpdate::Pointer);
itkSetMacro(r23Meta, R23MetaInformation::Pointer);
itkGetMacro(r23Meta, R23MetaInformation::Pointer);
itkSetMacro(Volume, InternalImageType::Pointer);
itkGetMacro(Volume, InternalImageType::Pointer);
itkSetMacro(PA, InternalImageType::Pointer);
itkGetMacro(PA, InternalImageType::Pointer);
itkSetMacro(LAT, InternalImageType::Pointer);
itkGetMacro(LAT, InternalImageType::Pointer);
itkSetMacro(InternalTransf1, InternalTransformMetaInformation::Pointer);
itkGetMacro(InternalTransf1, InternalTransformMetaInformation::Pointer);
itkSetMacro(InternalTransf2, InternalTransformMetaInformation::Pointer);
itkGetMacro(InternalTransf2, InternalTransformMetaInformation::Pointer);
itkSetMacro(filter1, ChangeInformationFilterType::Pointer);
itkGetMacro(filter1, ChangeInformationFilterType::Pointer);
itkSetMacro(filter2, ChangeInformationFilterType::Pointer);
itkGetMacro(filter2, ChangeInformationFilterType::Pointer);
itkSetMacro(interpolator1, InterpolatorType::Pointer);
itkGetMacro(interpolator1, InterpolatorType::Pointer);
itkSetMacro(interpolator2, InterpolatorType::Pointer);
itkGetMacro(interpolator2, InterpolatorType::Pointer);
itkSetMacro(TransformMetaInfo, TransformMetaInformation::Pointer);
itkGetMacro(TransformMetaInfo, TransformMetaInformation::Pointer);
/* ---- User provided END */
itkSetMacro(roiAutoReg1, RoiForRegistration);
itkGetMacro(roiAutoReg1, RoiForRegistration);
itkSetMacro(roiAutoReg2, RoiForRegistration);
itkGetMacro(roiAutoReg2, RoiForRegistration);
/** Auto Reg23 methods */
/** Initialize the registration pipeline*/
void InitializeRegistration(double stepLength,
double maxTranslation,
eDegreeOfFreedomType dof);
/** Start the registration process*/
int StartRegistration(std::string extraInfo);
/** Get the current cost function value from the optimizer*/
double GetOptimizerValue();
/** Get the cost function value for the current transform*/
double GetCurrentMetricValue();
/** Get the cost function value for the current transform*/
Optimizer::ParametersType GetCurrentPosition();
/** Auto Reg23 methods END */
protected:
itkReg23();
virtual ~itkReg23();
void PrintSelf(std::ostream& os, itk::Indent indent) const;
private:
itkReg23(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
/** Optimizer which tries to find the minimun (Powell Optimizer)*/
using PowellOptimizerType = itk::PowellOptimizer;
/** 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>;
RegistrationType::Pointer
registration;
MetricType::Pointer
NCCmetric;
MIMetricType::Pointer
MImetric;
PowellOptimizerType::Pointer
PowellOptimizer;
AmoebaOptimizerType::Pointer
AmoebaOptimizer;
ExhaustiveOptimizerType::Pointer
ExhaustiveOptimizer;
/* ---- User provided */
R23MetaInformation::Pointer
m_r23Meta;
InternalImageType::Pointer
m_Volume,
m_PA,
m_LAT;
InternalTransformMetaInformation::Pointer
m_InternalTransf1,
m_InternalTransf2;
ChangeInformationFilterType::Pointer
m_filter1,
m_filter2;
InterpolatorType::Pointer
m_interpolator1,
m_interpolator2;
TransformMetaInformation::Pointer
m_TransformMetaInfo;
CommandIterationUpdate::Pointer
m_OptimizerObserver;
/* ---- User provided END */
TransformType::Pointer
IsocTransform,
Transform1,
Transform2;
RoiForRegistration
m_roiAutoReg1,
m_roiAutoReg2;
bool m_UseDumptoFile;
double m_OptmizerValue;
};
}
#endif

5
reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.h
View File

@ -182,6 +182,10 @@ public:
itkSetMacro(Threshold, double);
itkGetMacro(Threshold, double);
/** Set and get the Panel Offset */
itkSetMacro(PanelOffset, double);
itkGetMacro(PanelOffset, double);
/** Check if a point is inside the image buffer.
* \warning For efficiency, no validity checking of
* the input image pointer is done. */
@ -232,6 +236,7 @@ protected:
double m_Threshold;
double m_FocalPointToIsocenterDistance; // Focal point to isocenter distance
double m_ProjectionAngle; // Linac gantry rotation angle in radians
double m_PanelOffset;
private:
void

20
reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx
View File

@ -62,6 +62,7 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::gSiddonJac
m_FocalPointToIsocenterDistance = 1000.; // Focal point to isocenter distance in mm.
m_ProjectionAngle = 0.; // Angle in radians betweeen projection central axis and reference axis
m_Threshold = 0.; // Intensity threshold, below which is ignored.
m_PanelOffset = 0.;
m_SourcePoint[0] = 0.;
m_SourcePoint[1] = 0.;
@ -148,14 +149,15 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
// m_SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
}
PointType PointReq = point;
PointReq[0] += m_PanelOffset;
drrPixelWorld = m_InverseTransform->TransformPoint(point);
drrPixelWorld = m_InverseTransform->TransformPoint(PointReq);
PointType SlidingSourcePoint = m_SourcePoint;
SlidingSourcePoint[0] = 0.;
SlidingSourcePoint[1] = point[1];
SlidingSourcePoint[2] = 0.;
//PointType SourceWorld = m_InverseTransform->TransformPoint(m_SourcePoint);
PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);
//std::cout<<"SourceWorld: "<<SourceWorld<<std::endl;
@ -173,8 +175,22 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
regionCT = inputPtr->GetLargestPossibleRegion();
sizeCT = regionCT.GetSize();
// If Pixel position (in mm) is outside bounds of CT (zero-based)
// assign 0 at the pixel and move on. Ensures regular spacing of resulting
// DRR
float xSizeCT = sizeCT[0] * ctPixelSpacing[0];
float ySizeCT = sizeCT[1] * ctPixelSpacing[1];
float zSizeCT = sizeCT[2] * ctPixelSpacing[2];
float xDrrPix = drrPixelWorld[0];float yDrrPix = drrPixelWorld[1];float zDrrPix = drrPixelWorld[2];
// if(zDrrPix < 0 /*|| yDrrPix < 0 || xDrrPix < 0*/)
// std::cout << drrPixelWorld[0]<<" " <<drrPixelWorld[1]<<" "<<drrPixelWorld[2]<<" / ";
if(xDrrPix>xSizeCT || yDrrPix>ySizeCT || zDrrPix>zSizeCT){
pixval = static_cast<OutputType>(0.0);
return pixval;
}
// calculate the detector position for this pixel center by moving
// 2*m_FocalPointToIsocenterDistance from the source in the pixel
// directions

19
reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.cxx
View File

@ -66,6 +66,22 @@ vtkContourTopogramProjectionFilter::~vtkContourTopogramProjectionFilter()
}
}
void vtkContourTopogramProjectionFilter::cleanUpFilter(){
this->dProjectionLPSOffset[0]=0.;
this->dProjectionLPSOffset[1]=0.;
this->dProjectionLPSOffset[2]=0.;
this->dImportOffsetLPS[0] = 0.;
this->dImportOffsetLPS[1] = 0.;
this->dImportOffsetLPS[2] = 0.;
this->dSCD = 0.;
this->m_RefTransform = nullptr;
this->m_Transform = nullptr;
}
int vtkContourTopogramProjectionFilter::RequestData(
vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
@ -94,6 +110,7 @@ int vtkContourTopogramProjectionFilter::RequestData(
// If no points, then nothing to do.
if (points == nullptr)
{
std::cout << "Cannot Project; no input points" << std::endl;
vtkDebugMacro("Cannot Project; no input points");
return 1;
}
@ -101,6 +118,7 @@ int vtkContourTopogramProjectionFilter::RequestData(
// If reference transform, then nothing to do.
if (m_RefTransform == nullptr)
{
std::cout << "Cannot Project; no input reference projection transform" << std::endl;
vtkDebugMacro("Cannot Project; no input reference projection transform");
return 1;
}
@ -108,6 +126,7 @@ int vtkContourTopogramProjectionFilter::RequestData(
// If transform, then nothing to do.
if (m_Transform == nullptr)
{
std::cout << "Cannot Project; no input projection transform" << std::endl;
vtkDebugMacro("Cannot Project; no input projection transform");
return 1;
}

2
reg23Topograms/itkDTRrecon/vtkContourTopogramProjectionFilter.h
View File

@ -41,6 +41,8 @@ public:
void SetSCD(const double dVal);
void SetImportOffsetLPS(const double *);
void cleanUpFilter();
protected:
vtkContourTopogramProjectionFilter();
~vtkContourTopogramProjectionFilter() override;