Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions documentation/docs/rtk_3_migration_guide.md
Original file line number Diff line number Diff line change
Expand Up @@ -43,3 +43,8 @@ Most users intuitively expect this for the backprojector matched to `rtk::CudaFo
## CUDA forward / ray-cast step size default

The CUDA forward projector and the CUDA ray-cast backprojector now use the minimum voxel spacing of the input volume as the ray step when StepSize is not set. The previous default value was 1.

## Positions now use `itk::Point` instead of `itk::Vector`

APIs that represent physical *positions* were changed from `itk::Vector<double,3>` to `itk::Point<double,3>`. This makes a clear distinction between points (locations) and vectors (directions or offsets).
Projection iterators (`rtkProjectionsRegionConstIteratorRayBased*`), Joseph forward/back projectors, ray-box filters, convex/quadric shapes, Ora geometry readers, Forbild phantom readers, and test/phantom helpers now expect `PointType` for positions (source, pixel, box corners, shape centers).
2 changes: 1 addition & 1 deletion examples/FirstReconstruction/FirstCudaReconstruction.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@ main(int argc, char ** argv)
REIType::Pointer rei = REIType::New();
rei->SetDensity(2.);
rei->SetAngle(0.);
rei->SetCenter(itk::MakeVector(0., 0., 10.));
rei->SetCenter(itk::MakePoint(0., 0., 10.));
rei->SetAxis(itk::MakeVector(50., 50., 50.));
rei->SetGeometry(geometry);
rei->SetInput(constantImageSource->GetOutput());
Expand Down
2 changes: 1 addition & 1 deletion examples/FirstReconstruction/FirstReconstruction.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ main(int argc, char ** argv)
REIType::Pointer rei = REIType::New();
rei->SetDensity(2.);
rei->SetAngle(0.);
rei->SetCenter(itk::MakeVector(0., 0., 10.));
rei->SetCenter(itk::MakePoint(0., 0., 10.));
rei->SetAxis(itk::MakeVector(50., 50., 50.));
rei->SetGeometry(geometry);
rei->SetInput(constantImageSource->GetOutput());
Expand Down
2 changes: 1 addition & 1 deletion include/rtkConvexShape.h
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ class RTK_EXPORT ConvexShape : public itk::DataObject
/** Convenient type alias. */
static constexpr unsigned int Dimension = 3;
using ScalarType = double;
using PointType = itk::Vector<ScalarType, Dimension>;
using PointType = itk::Point<ScalarType, Dimension>;
using VectorType = itk::Vector<ScalarType, Dimension>;
using RotationMatrixType = itk::Matrix<ScalarType, Dimension, Dimension>;

Expand Down
12 changes: 6 additions & 6 deletions include/rtkDrawBoxImageFilter.h
Original file line number Diff line number Diff line change
Expand Up @@ -77,10 +77,10 @@ class ITK_TEMPLATE_EXPORT DrawBoxImageFilter : public DrawConvexImageFilter<TInp
SetBoxFromImage(const ImageBaseType * img, bool bWithExternalHalfPixelBorder = true);

/** Get/Set the box parameters. See rtk::BoxShape. */
itkGetMacro(BoxMin, VectorType);
itkSetMacro(BoxMin, VectorType);
itkGetMacro(BoxMax, VectorType);
itkSetMacro(BoxMax, VectorType);
itkGetMacro(BoxMin, PointType);
itkSetMacro(BoxMin, PointType);
itkGetMacro(BoxMax, PointType);
itkSetMacro(BoxMax, PointType);
itkGetMacro(Direction, RotationMatrixType);
itkSetMacro(Direction, RotationMatrixType);

Expand All @@ -96,8 +96,8 @@ class ITK_TEMPLATE_EXPORT DrawBoxImageFilter : public DrawConvexImageFilter<TInp
std::vector<VectorType> m_PlaneDirections;
std::vector<ScalarType> m_PlanePositions;

VectorType m_BoxMin{ 0. };
VectorType m_BoxMax{ 0. };
PointType m_BoxMin{ 0. };
PointType m_BoxMax{ 0. };
RotationMatrixType m_Direction;
};

Expand Down
24 changes: 10 additions & 14 deletions include/rtkFourDSARTConeBeamReconstructionFilter.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -236,20 +236,16 @@ FourDSARTConeBeamReconstructionFilter<VolumeSeriesType, ProjectionStackType>::Ge

// Create the m_RayBoxFiltersectionImageFilter
m_RayBoxFilter->SetGeometry(this->GetGeometry());
itk::Vector<double, 3> Corner1, Corner2;

Corner1[0] = this->GetInput(0)->GetOrigin()[0];
Corner1[1] = this->GetInput(0)->GetOrigin()[1];
Corner1[2] = this->GetInput(0)->GetOrigin()[2];
Corner2[0] = this->GetInput(0)->GetOrigin()[0] +
this->GetInput(0)->GetLargestPossibleRegion().GetSize()[0] * this->GetInput(0)->GetSpacing()[0];
Corner2[1] = this->GetInput(0)->GetOrigin()[1] +
this->GetInput(0)->GetLargestPossibleRegion().GetSize()[1] * this->GetInput(0)->GetSpacing()[1];
Corner2[2] = this->GetInput(0)->GetOrigin()[2] +
this->GetInput(0)->GetLargestPossibleRegion().GetSize()[2] * this->GetInput(0)->GetSpacing()[2];

m_RayBoxFilter->SetBoxMin(Corner1);
m_RayBoxFilter->SetBoxMax(Corner2);
itk::Point<double, 3> corner1, corner2;
for (unsigned int i = 0; i < 3; ++i)
{
corner1[i] = this->GetInput(0)->GetOrigin()[i];
corner2[i] = this->GetInput(0)->GetOrigin()[i] +
this->GetInput(0)->GetLargestPossibleRegion().GetSize()[i] * this->GetInput(0)->GetSpacing()[i];
}

m_RayBoxFilter->SetBoxMin(corner1);
m_RayBoxFilter->SetBoxMax(corner2);

m_RayBoxFilter->UpdateOutputInformation();
m_ExtractFilter->UpdateOutputInformation();
Expand Down
22 changes: 10 additions & 12 deletions include/rtkIterativeFDKConeBeamReconstructionFilter.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -97,18 +97,16 @@ IterativeFDKConeBeamReconstructionFilter<TInputImage, TOutputImage, TFFTPrecisio
m_ConstantProjectionStackSource->SetConstant(0);

// Set box for the m_RayBoxFiltersectionImageFilter
itk::Vector<double, 3> Corner1, Corner2;
Corner1[0] = this->GetInput(0)->GetOrigin()[0];
Corner1[1] = this->GetInput(0)->GetOrigin()[1];
Corner1[2] = this->GetInput(0)->GetOrigin()[2];
Corner2[0] = this->GetInput(0)->GetOrigin()[0] +
this->GetInput(0)->GetLargestPossibleRegion().GetSize()[0] * this->GetInput(0)->GetSpacing()[0];
Corner2[1] = this->GetInput(0)->GetOrigin()[1] +
this->GetInput(0)->GetLargestPossibleRegion().GetSize()[1] * this->GetInput(0)->GetSpacing()[1];
Corner2[2] = this->GetInput(0)->GetOrigin()[2] +
this->GetInput(0)->GetLargestPossibleRegion().GetSize()[2] * this->GetInput(0)->GetSpacing()[2];
m_RayBoxFilter->SetBoxMin(Corner1);
m_RayBoxFilter->SetBoxMax(Corner2);
itk::Point<double, 3> corner1, corner2;
for (unsigned int i = 0; i < 3; ++i)
{
corner1[i] = this->GetInput(0)->GetOrigin()[i];
corner2[i] = this->GetInput(0)->GetOrigin()[i] +
this->GetInput(0)->GetLargestPossibleRegion().GetSize()[i] * this->GetInput(0)->GetSpacing()[i];
}

m_RayBoxFilter->SetBoxMin(corner1);
m_RayBoxFilter->SetBoxMax(corner2);

// Initial internal connections
m_DisplacedDetectorFilter->SetInput(this->GetInput(1));
Expand Down
11 changes: 6 additions & 5 deletions include/rtkJosephBackProjectionImageFilter.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -109,8 +109,8 @@ JosephBackProjectionImageFilter<TInputImage,
itIn = InputRegionIterator::New(this->GetInput(1), buffReg, geometry, volPPToIndex);

// Create intersection functions, one for each possible main direction
auto box = BoxShape::New();
typename BoxShape::VectorType boxMin, boxMax;
auto box = BoxShape::New();
typename BoxShape::PointType boxMin, boxMax;
for (unsigned int i = 0; i < Dimension; i++)
{
boxMin[i] = this->GetOutput()->GetRequestedRegion().GetIndex()[i];
Expand All @@ -127,11 +127,12 @@ JosephBackProjectionImageFilter<TInputImage,
double superiorClip = 1. - m_InferiorClip;

// Go over each pixel of the projection
typename BoxShape::VectorType stepMM, np, fp;
typename BoxShape::VectorType stepMM;
typename BoxShape::PointType np, fp;
for (unsigned int pix = 0; pix < buffReg.GetNumberOfPixels(); pix++, itIn->Next())
{
typename InputRegionIterator::PointType pixelPosition = itIn->GetPixelPosition();
typename InputRegionIterator::PointType dirVox = -itIn->GetSourceToPixel();
typename InputRegionIterator::PointType pixelPosition = itIn->GetPixelPosition();
typename InputRegionIterator::VectorType dirVox = -itIn->GetSourceToPixel();

// Select main direction
unsigned int mainDir = 0;
Expand Down
7 changes: 4 additions & 3 deletions include/rtkJosephForwardAttenuatedProjectionImageFilter.h
Original file line number Diff line number Diff line change
Expand Up @@ -189,6 +189,7 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttenuated
{
public:
using VectorType = itk::Vector<double, 3>;
using PointType = itk::Point<double, 3>;

ProjectedValueAccumulationAttenuated() = default;
~ProjectedValueAccumulationAttenuated() = default;
Expand All @@ -210,10 +211,10 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulationAttenuated
TOutput & output,
const TOutput & rayCastValue,
const VectorType & /*stepInMM*/,
const VectorType & itkNotUsed(source),
const PointType & itkNotUsed(source),
const VectorType & itkNotUsed(sourceToPixel),
const VectorType & itkNotUsed(nearestPoint),
const VectorType & itkNotUsed(farthestPoint))
const PointType & itkNotUsed(nearestPoint),
const PointType & itkNotUsed(farthestPoint))
{
output = input + rayCastValue;
m_Attenuation[threadId] = 0;
Expand Down
7 changes: 4 additions & 3 deletions include/rtkJosephForwardProjectionImageFilter.h
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,7 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulation
{
public:
using VectorType = itk::Vector<double, 3>;
using PointType = itk::Point<double, 3>;

ProjectedValueAccumulation() = default;
~ProjectedValueAccumulation() = default;
Expand All @@ -135,10 +136,10 @@ class ITK_TEMPLATE_EXPORT ProjectedValueAccumulation
TOutput & output,
const TOutput & rayCastValue,
const VectorType & stepInMM,
const VectorType & itkNotUsed(source),
const PointType & itkNotUsed(source),
const VectorType & itkNotUsed(sourceToPixel),
const VectorType & itkNotUsed(nearestPoint),
const VectorType & itkNotUsed(farthestPoint)) const
const PointType & itkNotUsed(nearestPoint),
const PointType & itkNotUsed(farthestPoint)) const
{
output = input + rayCastValue * stepInMM.GetNorm();
}
Expand Down
14 changes: 8 additions & 6 deletions include/rtkJosephForwardProjectionImageFilter.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -221,8 +221,8 @@ JosephForwardProjectionImageFilter<TInputImage,
}

// Create intersection functions, one for each possible main direction
auto box = BoxShape::New();
typename BoxShape::VectorType boxMin, boxMax;
auto box = BoxShape::New();
typename BoxShape::PointType boxMin, boxMax;
for (unsigned int i = 0; i < Dimension; i++)
{
boxMin[i] = this->GetInput(1)->GetBufferedRegion().GetIndex()[i];
Expand All @@ -239,11 +239,13 @@ JosephForwardProjectionImageFilter<TInputImage,
double superiorClip = 1. - m_InferiorClip;

// Go over each pixel of the projection
typename BoxShape::VectorType stepMM, np, fp;
typename BoxShape::VectorType stepMM;
typename InputRegionIterator::PointType np, fp;
for (unsigned int pix = 0; pix < outputRegionForThread.GetNumberOfPixels(); pix++, itIn->Next(), ++itOut)
{
typename InputRegionIterator::PointType pixelPosition = itIn->GetPixelPosition();
typename InputRegionIterator::PointType dirVox = -itIn->GetSourceToPixel();
stepMM.Fill(0.);
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems to fix the failing test that occurred because of pixelPosition that was corrected to stepMM line 422

typename InputRegionIterator::PointType pixelPosition = itIn->GetPixelPosition();
typename InputRegionIterator::VectorType dirVox = -itIn->GetSourceToPixel();

// Select main direction
unsigned int mainDir = 0;
Expand Down Expand Up @@ -417,7 +419,7 @@ JosephForwardProjectionImageFilter<TInputImage,
}
else
m_ProjectedValueAccumulation(
threadId, itIn->Get(), itOut.Value(), {}, pixelPosition, pixelPosition, dirVox, pixelPosition, pixelPosition);
threadId, itIn->Get(), itOut.Value(), {}, stepMM, pixelPosition, dirVox, pixelPosition, pixelPosition);
}
delete itIn;
}
Expand Down
18 changes: 10 additions & 8 deletions include/rtkMaskCollimationImageFilter.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -73,26 +73,28 @@ MaskCollimationImageFilter<TInputImage, TOutputImage>::DynamicThreadedGenerateDa
for (unsigned int pix = 0; pix < npixperslice; pix++, itIn->Next(), ++itOut)
{
using PointType = typename InputRegionIterator::PointType;
PointType source = itIn->GetSourcePosition();
double sourceNorm = source.GetNorm();
PointType pixel = itIn->GetPixelPosition();
using VectorType = typename InputRegionIterator::VectorType;

PointType pixelToSource(source - pixel);
PointType sourceDir = source / sourceNorm;
VectorType source = itIn->GetSourcePosition().GetVectorFromOrigin();
double sourceNorm = source.GetNorm();
VectorType pixel = itIn->GetPixelPosition().GetVectorFromOrigin();

VectorType pixelToSource(source - pixel);
VectorType sourceDir = source / sourceNorm;

// Get the 3D position of the ray and the main plane
// (at 0,0,0 and orthogonal to source to center)
const double mag = sourceNorm / (pixelToSource * sourceDir);
PointType intersection = source - mag * pixelToSource;
VectorType intersection = source - mag * pixelToSource;

// The first coordinate on the plane is measured along the detector
// orthogonal to the plane defined by the detector source direction
// and the secon detector vector
PointType v;
VectorType v;
v[0] = this->GetGeometry()->GetRotationMatrix(iProj)[m_RotationAxisIndex][0];
v[1] = this->GetGeometry()->GetRotationMatrix(iProj)[m_RotationAxisIndex][1];
v[2] = this->GetGeometry()->GetRotationMatrix(iProj)[m_RotationAxisIndex][2];
PointType u = CrossProduct(v, sourceDir);
VectorType u = CrossProduct(v, sourceDir);
u /= u.GetNorm();
double ucoord = u * intersection;

Expand Down
7 changes: 4 additions & 3 deletions include/rtkMaximumIntensityProjectionImageFilter.h
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,7 @@ class ITK_TEMPLATE_EXPORT MaximumIntensityProjectedValueAccumulation
{
public:
using VectorType = itk::Vector<double, 3>;
using PointType = itk::Point<double, 3>;

bool
operator!=(const MaximumIntensityProjectedValueAccumulation &) const
Expand All @@ -96,10 +97,10 @@ class ITK_TEMPLATE_EXPORT MaximumIntensityProjectedValueAccumulation
TOutput & output,
const TOutput & rayCastValue,
const VectorType & stepInMM,
const VectorType & itkNotUsed(source),
const PointType & itkNotUsed(source),
const VectorType & itkNotUsed(sourceToPixel),
const VectorType & itkNotUsed(nearestPoint),
const VectorType & itkNotUsed(farthestPoint)) const
const PointType & itkNotUsed(nearestPoint),
const PointType & itkNotUsed(farthestPoint)) const
{
TOutput tmp = static_cast<TOutput>(input);
if (tmp < rayCastValue)
Expand Down
9 changes: 5 additions & 4 deletions include/rtkProjectionsRegionConstIteratorRayBased.h
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,8 @@ class ITK_TEMPLATE_EXPORT ProjectionsRegionConstIteratorRayBased : public itk::I
/** Types inherited from the Superclass */
using OffsetValueType = typename Superclass::OffsetValueType;
using RegionType = typename Superclass::RegionType;
using PointType = typename itk::Vector<double, 3>;
using VectorType = typename itk::Vector<double, 3>;
using PointType = typename itk::Point<double, 3>;
using IndexValueType = typename Superclass::IndexValueType;

using MatrixType = itk::Matrix<double, 3, 4>;
Expand Down Expand Up @@ -130,15 +131,15 @@ class ITK_TEMPLATE_EXPORT ProjectionsRegionConstIteratorRayBased : public itk::I
{
return this->m_PixelPosition;
}
const PointType &
const VectorType &
GetSourceToPixel()
{
return this->m_SourceToPixel;
}

/** Computes and returns a unit vector pointing from the source to the
* current pixel, i.e., GetSourceToPixel()/||GetSourceToPixel()||. */
const PointType
const VectorType
GetDirection()
{
return m_SourceToPixel / m_SourceToPixel.GetNorm();
Expand All @@ -158,7 +159,7 @@ class ITK_TEMPLATE_EXPORT ProjectionsRegionConstIteratorRayBased : public itk::I
MatrixType m_PostMultiplyMatrix;
PointType m_SourcePosition;
PointType m_PixelPosition;
PointType m_SourceToPixel;
VectorType m_SourceToPixel;
};
} // namespace rtk

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,8 @@ class ITK_TEMPLATE_EXPORT ProjectionsRegionConstIteratorRayBasedWithCylindricalP
using MatrixType = typename Superclass::MatrixType;
using IndexValueType = typename Superclass::IndexValueType;

using PointType = typename itk::Vector<double, 3>;
using PointType = typename itk::Point<double, 3>;
using VectorType = typename itk::Vector<double, 3>;
using HomogeneousMatrixType = itk::Matrix<double, 4, 4>;

/** Constructor establishes an iterator to walk a particular image and a
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,9 @@ ProjectionsRegionConstIteratorRayBasedWithCylindricalPanel<TImage>::NewProjectio
// Set source position in volume indices
// GetSourcePosition() returns coordinates in mm. Multiplying by
// volPPToIndex gives the corresponding volume index
this->m_SourcePosition = this->m_PostMultiplyMatrix * this->m_Geometry->GetSourcePosition(iProj);
VectorType sourcePos = this->m_PostMultiplyMatrix * this->m_Geometry->GetSourcePosition(iProj);
for (unsigned int i = 0; i < this->GetImageDimension(); ++i)
this->m_SourcePosition[i] = sourcePos[i];

// Compute matrix to transform projection index to position on a flat panel
// if the panel were flat, before accounting for the curvature
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,8 @@ class ITK_TEMPLATE_EXPORT ProjectionsRegionConstIteratorRayBasedWithFlatPanel
/** Types inherited from the Superclass */
using OffsetValueType = typename Superclass::OffsetValueType;
using RegionType = typename Superclass::RegionType;
using PointType = typename itk::Vector<double, 3>;
using PointType = typename itk::Point<double, 3>;
using VectorType = typename itk::Vector<double, 3>;
using MatrixType = typename Superclass::MatrixType;
using HomogeneousMatrixType = itk::Matrix<double, 4, 4>;

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,9 @@ ProjectionsRegionConstIteratorRayBasedWithFlatPanel<TImage>::NewProjection()
// Set source position in volume indices
// GetSourcePosition() returns coordinates in mm. Multiplying by
// volPPToIndex gives the corresponding volume index
this->m_SourcePosition = this->m_PostMultiplyMatrix * this->m_Geometry->GetSourcePosition(this->m_PositionIndex[2]);
VectorType transformed = this->m_PostMultiplyMatrix * this->m_Geometry->GetSourcePosition(this->m_PositionIndex[2]);
for (unsigned int i = 0; i < this->GetImageDimension(); ++i)
this->m_SourcePosition[i] = transformed[i];

// Compute matrix to transform projection index to volume index
// IndexToPhysicalPointMatrix maps the 2D index of a projection's pixel to its 2D position on the detector (in mm)
Expand Down
10 changes: 6 additions & 4 deletions src/rtkBoxShape.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -145,17 +145,19 @@ BoxShape ::SetBoxFromImage(const ImageBaseType * img, bool bWithExternalHalfPixe
itkGenericExceptionMacro(<< "BoxShape and image dimensions must agree");

// BoxShape corner 1
m_BoxMin = img->GetOrigin().GetVectorFromOrigin();
m_BoxMin = img->GetOrigin();
if (bWithExternalHalfPixelBorder)
m_BoxMin -= img->GetDirection() * img->GetSpacing() * 0.5;

// BoxShape corner 2
VectorType max;
for (unsigned int i = 0; i < Dimension; i++)
if (bWithExternalHalfPixelBorder)
m_BoxMax[i] = img->GetSpacing()[i] * img->GetLargestPossibleRegion().GetSize()[i];
max[i] = img->GetSpacing()[i] * img->GetLargestPossibleRegion().GetSize()[i];
else
m_BoxMax[i] = img->GetSpacing()[i] * (img->GetLargestPossibleRegion().GetSize()[i] - 1);
m_BoxMax = m_BoxMin + img->GetDirection() * m_BoxMax;
max[i] = img->GetSpacing()[i] * (img->GetLargestPossibleRegion().GetSize()[i] - 1);
max = img->GetDirection() * max;
m_BoxMax = m_BoxMin + max;

SetDirection(img->GetDirection());
}
Expand Down
Loading
Loading