From eeb2c919a58064a4811efdccc9703f55cf1f133c Mon Sep 17 00:00:00 2001 From: Cyril Mory Date: Thu, 8 Jan 2026 11:16:46 +0100 Subject: [PATCH 1/3] ENH: Move 3D+time test data generation code into a helper file Move the 3D+time test data generation from test/rtkfourdroostertest.cxx into test/rtkFourDTestHelper.h and .hxx Modify test/rtkfourdroostertest.cxx accordingly --- test/rtkFourDTestHelper.h | 78 ++++++++++ test/rtkFourDTestHelper.hxx | 277 +++++++++++++++++++++++++++++++++++ test/rtkfourdroostertest.cxx | 264 +++------------------------------ 3 files changed, 372 insertions(+), 247 deletions(-) create mode 100644 test/rtkFourDTestHelper.h create mode 100644 test/rtkFourDTestHelper.hxx diff --git a/test/rtkFourDTestHelper.h b/test/rtkFourDTestHelper.h new file mode 100644 index 000000000..bbae82b82 --- /dev/null +++ b/test/rtkFourDTestHelper.h @@ -0,0 +1,78 @@ +/*========================================================================= + * + * Copyright RTK Consortium + * + * 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 + * + * https://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 rtkFourDTestHelper_h +#define rtkFourDTestHelper_h + +#include + +#include +#include + +#ifdef RTK_USE_CUDA +# include +#endif + +#include "rtkThreeDCircularProjectionGeometry.h" + +namespace rtk +{ + +/** Structure containing all the test data required for 3D+time filters */ +template +struct FourDTestData +{ + using DVFVectorType = itk::CovariantVector; + +#ifdef RTK_USE_CUDA + using VolumeSeriesType = itk::CudaImage; + using ProjectionStackType = itk::CudaImage; + using VolumeType = itk::CudaImage; + using DVFSequenceImageType = itk::CudaImage; +#else + using VolumeSeriesType = itk::Image; + using ProjectionStackType = itk::Image; + using VolumeType = itk::Image; + using DVFSequenceImageType = itk::Image; +#endif + + typename VolumeType::Pointer SingleVolume; + typename VolumeSeriesType::Pointer InitialVolumeSeries; + typename ProjectionStackType::Pointer Projections; + typename rtk::ThreeDCircularProjectionGeometry::Pointer Geometry; + + typename DVFSequenceImageType::Pointer DVF; + typename DVFSequenceImageType::Pointer InverseDVF; + + typename VolumeSeriesType::Pointer GroundTruth; + + std::string SignalFileName; + std::vector Signal; +}; + +/** Generates a full 3D+time dataset (phantom, projections, signal, DVF, ground truth) */ +template +FourDTestData +GenerateFourDTestData(bool fastTests); + +} // namespace rtk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "rtkFourDTestHelper.hxx" +#endif + +#endif diff --git a/test/rtkFourDTestHelper.hxx b/test/rtkFourDTestHelper.hxx new file mode 100644 index 000000000..475335e80 --- /dev/null +++ b/test/rtkFourDTestHelper.hxx @@ -0,0 +1,277 @@ +/*========================================================================= + * + * Copyright RTK Consortium + * + * 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 + * + * https://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 rtkFourDTestHelper_hxx +#define rtkFourDTestHelper_hxx + +#include + +#include +#include +#include +#include + +#include "rtkConstantImageSource.h" +#include "rtkRayEllipsoidIntersectionImageFilter.h" +#include "rtkDrawEllipsoidImageFilter.h" +#include "rtkTest.h" + +namespace rtk +{ + +template +FourDTestData +GenerateFourDTestData(bool fastTests) +{ + using DVFVectorType = itk::CovariantVector; + +#ifdef RTK_USE_CUDA + using VolumeSeriesType = itk::CudaImage; + using ProjectionStackType = itk::CudaImage; + using VolumeType = itk::CudaImage; + using DVFSequenceImageType = itk::CudaImage; +#else + using VolumeSeriesType = itk::Image; + using ProjectionStackType = itk::Image; + using VolumeType = itk::Image; + using DVFSequenceImageType = itk::Image; +#endif + using GeometryType = rtk::ThreeDCircularProjectionGeometry; + +#if FAST_TESTS_NO_CHECKS + constexpr unsigned int NumberOfProjectionImages = 5; +#else + constexpr unsigned int NumberOfProjectionImages = 64; +#endif + + FourDTestData data; + + /* ========================= + * Constant volume sources + * ========================= */ + using ConstantVolumeSourceType = rtk::ConstantImageSource; + auto volumeSource = ConstantVolumeSourceType::New(); + + auto origin = itk::MakePoint(-63., -31., -63.); + auto spacing = fastTests ? itk::MakeVector(16., 8., 16.) : itk::MakeVector(4., 4., 4.); + auto size = fastTests ? itk::MakeSize(8, 8, 8) : itk::MakeSize(32, 16, 32); + + volumeSource->SetOrigin(origin); + volumeSource->SetSpacing(spacing); + volumeSource->SetSize(size); + volumeSource->SetConstant(0.); + volumeSource->Update(); + + data.SingleVolume = volumeSource->GetOutput(); + + using ConstantFourDSourceType = rtk::ConstantImageSource; + auto fourDSource = ConstantFourDSourceType::New(); + + auto fourDOrigin = itk::MakePoint(-63., -31., -63., 0.); + auto fourDSpacing = fastTests ? itk::MakeVector(16., 8., 16., 1.) : itk::MakeVector(4., 4., 4., 1.); + auto fourDSize = fastTests ? itk::MakeSize(8, 8, 8, 2) : itk::MakeSize(32, 16, 32, 8); + + fourDSource->SetOrigin(fourDOrigin); + fourDSource->SetSpacing(fourDSpacing); + fourDSource->SetSize(fourDSize); + fourDSource->SetConstant(0.); + fourDSource->Update(); + + data.InitialVolumeSeries = fourDSource->GetOutput(); + + /* ========================= + * Projection accumulation + * ========================= */ + using ConstantProjectionSourceType = rtk::ConstantImageSource; + auto projectionsSource = ConstantProjectionSourceType::New(); + + auto projOrigin = itk::MakePoint(-254., -254., -254.); + auto projSpacing = fastTests ? itk::MakeVector(32., 32., 32.) : itk::MakeVector(8., 8., 1.); + auto projSize = + fastTests ? itk::MakeSize(32, 32, NumberOfProjectionImages) : itk::MakeSize(64, 64, NumberOfProjectionImages); + + projectionsSource->SetOrigin(projOrigin); + projectionsSource->SetSpacing(projSpacing); + projectionsSource->SetSize(projSize); + projectionsSource->SetConstant(0.); + projectionsSource->Update(); + + auto oneProjectionSource = ConstantProjectionSourceType::New(); + projSize[2] = 1; + oneProjectionSource->SetOrigin(projOrigin); + oneProjectionSource->SetSpacing(projSpacing); + oneProjectionSource->SetSize(projSize); + oneProjectionSource->SetConstant(0.); + + data.Geometry = GeometryType::New(); + + using REIType = rtk::RayEllipsoidIntersectionImageFilter; + using PasteType = itk::PasteImageFilter; + + /* Allocate explicit accumulator image */ + auto accumulated = ProjectionStackType::New(); + accumulated->SetOrigin(projectionsSource->GetOutput()->GetOrigin()); + accumulated->SetSpacing(projectionsSource->GetOutput()->GetSpacing()); + accumulated->SetDirection(projectionsSource->GetOutput()->GetDirection()); + accumulated->SetRegions(projectionsSource->GetOutput()->GetLargestPossibleRegion()); + accumulated->Allocate(); + accumulated->FillBuffer(0.); + + auto paste = PasteType::New(); + paste->SetDestinationImage(accumulated); + + itk::Index<3> destIndex; + destIndex.Fill(0); + + // Create the signal, as a file and as an std::vector +#ifdef USE_CUDA + data.SignalFileName = "signal_4D_cuda.txt"; +#else + data.SignalFileName = "signal_4D.txt"; +#endif + + data.Signal = std::vector(); + std::ofstream signalFile(data.SignalFileName.c_str()); + + for (unsigned int i = 0; i < NumberOfProjectionImages; ++i) + { + data.Geometry->AddProjection(600., 1200., i * 360. / NumberOfProjectionImages, 0, 0, 0, 0, 20, 15); + + auto geom = GeometryType::New(); + geom->AddProjection(600., 1200., i * 360. / NumberOfProjectionImages, 0, 0, 0, 0, 20, 15); + + auto e1 = REIType::New(); + e1->SetInput(oneProjectionSource->GetOutput()); + e1->SetGeometry(geom); + e1->SetDensity(2.); + e1->SetAxis(itk::MakeVector(60., 30., 60.)); + e1->SetCenter(itk::MakeVector(0., 0., 0.)); + e1->InPlaceOff(); + e1->Update(); + + auto e2 = REIType::New(); + e2->SetInput(e1->GetOutput()); + e2->SetGeometry(geom); + e2->SetDensity(-1.); + e2->SetAxis(itk::MakeVector(8., 8., 8.)); + auto center = itk::MakeVector(4 * (itk::Math::abs((4 + i) % 8 - 4.) - 2.), 0., 0.); + e2->SetCenter(center); + e2->InPlaceOff(); + e2->Update(); + + paste->SetSourceImage(e2->GetOutput()); + paste->SetSourceRegion(e2->GetOutput()->GetLargestPossibleRegion()); + paste->SetDestinationIndex(destIndex); + paste->Update(); + + accumulated = paste->GetOutput(); + paste->SetDestinationImage(accumulated); + + destIndex[2]++; + + data.Signal.push_back((i % 8) / 8.); + signalFile << data.Signal.back() << std::endl; + } + + signalFile.close(); + data.Projections = accumulated; + + /* ========================= + * DVF & inverse DVF + * ========================= */ + auto dvf = DVFSequenceImageType::New(); + auto idvf = DVFSequenceImageType::New(); + + typename DVFSequenceImageType::RegionType region; + auto dvfSize = itk::MakeSize(fourDSize[0], fourDSize[1], fourDSize[2], 2); + region.SetSize(dvfSize); + + dvf->SetRegions(region); + dvf->SetOrigin(fourDOrigin); + dvf->SetSpacing(fourDSpacing); + dvf->Allocate(); + + idvf->SetRegions(region); + idvf->SetOrigin(fourDOrigin); + idvf->SetSpacing(fourDSpacing); + idvf->Allocate(); + + itk::ImageRegionIteratorWithIndex it(dvf, region); + itk::ImageRegionIteratorWithIndex iit(idvf, region); + + DVFVectorType v; + typename DVFSequenceImageType::IndexType centerIndex; + centerIndex.Fill(0); + centerIndex[0] = dvfSize[0] / 2; + centerIndex[1] = dvfSize[1] / 2; + centerIndex[2] = dvfSize[2] / 2; + + for (; !it.IsAtEnd(); ++it, ++iit) + { + v.Fill(0.); + auto d = it.GetIndex() - centerIndex; + if (0.3 * d[0] * d[0] + d[1] * d[1] + d[2] * d[2] < 40) + v[0] = (it.GetIndex()[3] == 0 ? -8. : 8.); + it.Set(v); + iit.Set(-v); + } + + data.DVF = dvf; + data.InverseDVF = idvf; + + /* ========================= + * Ground truth + * ========================= */ + auto join = itk::JoinSeriesImageFilter::New(); + + for (unsigned int t = 0; t < fourDSize[3]; ++t) + { + using DEType = rtk::DrawEllipsoidImageFilter; + + auto de1 = DEType::New(); + de1->SetInput(volumeSource->GetOutput()); + de1->SetDensity(2.); + de1->SetAxis(itk::MakeVector(60., 30., 60.)); + de1->SetCenter(itk::MakeVector(0., 0., 0.)); + de1->InPlaceOff(); + de1->Update(); + + auto de2 = DEType::New(); + de2->SetInput(de1->GetOutput()); + de2->SetDensity(-1.); + de2->SetAxis(itk::MakeVector(8., 8., 8.)); + de2->SetCenter(itk::MakeVector(4 * (itk::Math::abs((4 + t) % 8 - 4.) - 2.), 0., 0.)); + de2->InPlaceOff(); + de2->Update(); + + using DuplicatorType = itk::ImageDuplicator; + auto duplicator = DuplicatorType::New(); + duplicator->SetInputImage(de2->GetOutput()); + duplicator->Update(); + + join->SetInput(t, duplicator->GetOutput()); + } + + join->Update(); + data.GroundTruth = join->GetOutput(); + + return data; +} + +} // namespace rtk + +#endif diff --git a/test/rtkfourdroostertest.cxx b/test/rtkfourdroostertest.cxx index b57c865b4..3e96007ed 100644 --- a/test/rtkfourdroostertest.cxx +++ b/test/rtkfourdroostertest.cxx @@ -1,15 +1,9 @@ -#include #include -#include #include "rtkTest.h" -#include "rtkRayEllipsoidIntersectionImageFilter.h" -#include "rtkDrawEllipsoidImageFilter.h" -#include "rtkConstantImageSource.h" -#include "rtkFieldOfViewImageFilter.h" -#include "rtkCyclicDeformationImageFilter.h" #include "rtkFourDROOSTERConeBeamReconstructionFilter.h" #include "rtkPhasesToInterpolationWeights.h" +#include "rtkFourDTestHelper.h" /** * \file rtkfourdroostertest.cxx @@ -29,7 +23,6 @@ int main(int, char **) { using OutputPixelType = float; - using DVFVectorType = itk::CovariantVector; #ifdef USE_CUDA @@ -44,233 +37,12 @@ main(int, char **) using DVFSequenceImageType = itk::Image; #endif -#if FAST_TESTS_NO_CHECKS - constexpr unsigned int NumberOfProjectionImages = 5; -#else - constexpr unsigned int NumberOfProjectionImages = 64; -#endif - - // Constant image sources - using ConstantImageSourceType = rtk::ConstantImageSource; - auto tomographySource = ConstantImageSourceType::New(); - auto origin = itk::MakePoint(-63., -31., -63.); -#if FAST_TESTS_NO_CHECKS - auto size = itk::MakeSize(8, 8, 8); - auto spacing = itk::MakeVector(16., 8., 16.); -#else - auto size = itk::MakeSize(32, 16, 32); - auto spacing = itk::MakeVector(4., 4., 4.); -#endif - tomographySource->SetOrigin(origin); - tomographySource->SetSpacing(spacing); - tomographySource->SetSize(size); - tomographySource->SetConstant(0.); - - auto fourdSource = rtk::ConstantImageSource::New(); - auto fourDOrigin = itk::MakePoint(-63., -31., -63., 0.); -#if FAST_TESTS_NO_CHECKS - auto fourDSize = itk::MakeSize(8, 8, 8, 2); - auto fourDSpacing = itk::MakeVector(16., 8., 16., 1.); -#else - auto fourDSize = itk::MakeSize(32, 16, 32, 8); - auto fourDSpacing = itk::MakeVector(4., 4., 4., 1.); -#endif - fourdSource->SetOrigin(fourDOrigin); - fourdSource->SetSpacing(fourDSpacing); - fourdSource->SetSize(fourDSize); - fourdSource->SetConstant(0.); - - auto projectionsSource = ConstantImageSourceType::New(); - origin = itk::MakePoint(-254., -254., -254.); -#if FAST_TESTS_NO_CHECKS - size = itk::MakeSize(32, 32, NumberOfProjectionImages); - spacing = itk::MakeVector(32., 32., 32.); -#else - size = itk::MakeSize(64, 64, NumberOfProjectionImages); - spacing = itk::MakeVector(8., 8., 1.); -#endif - projectionsSource->SetOrigin(origin); - projectionsSource->SetSpacing(spacing); - projectionsSource->SetSize(size); - projectionsSource->SetConstant(0.); - - auto oneProjectionSource = ConstantImageSourceType::New(); - size[2] = 1; - oneProjectionSource->SetOrigin(origin); - oneProjectionSource->SetSpacing(spacing); - oneProjectionSource->SetSize(size); - oneProjectionSource->SetConstant(0.); - - // Geometry object - using GeometryType = rtk::ThreeDCircularProjectionGeometry; - auto geometry = GeometryType::New(); - - // Projections - using REIType = rtk::RayEllipsoidIntersectionImageFilter; - auto destinationIndex = itk::MakeIndex(0, 0, 0); - auto pasteFilter = itk::PasteImageFilter::New(); - pasteFilter->SetDestinationImage(projectionsSource->GetOutput()); - -#ifdef USE_CUDA - std::string signalFileName = "signal_4DRooster_cuda.txt"; -#else - std::string signalFileName = "signal_4DRooster.txt"; -#endif - - std::ofstream signalFile(signalFileName.c_str()); - for (unsigned int noProj = 0; noProj < NumberOfProjectionImages; noProj++) - { - geometry->AddProjection(600., 1200., noProj * 360. / NumberOfProjectionImages, 0, 0, 0, 0, 20, 15); - - // Geometry object - auto oneProjGeometry = GeometryType::New(); - oneProjGeometry->AddProjection(600., 1200., noProj * 360. / NumberOfProjectionImages, 0, 0, 0, 0, 20, 15); - - // Ellipse 1 - auto e1 = REIType::New(); - auto semiprincipalaxis = itk::MakeVector(60., 30., 60.); - auto center = itk::MakeVector(0., 0., 0.); - e1->SetInput(oneProjectionSource->GetOutput()); - e1->SetGeometry(oneProjGeometry); - e1->SetDensity(2.); - e1->SetAxis(semiprincipalaxis); - e1->SetCenter(center); - e1->SetAngle(0.); - e1->InPlaceOff(); - e1->Update(); - - // Ellipse 2 - auto e2 = REIType::New(); - semiprincipalaxis.Fill(8.); - center[0] = 4 * (itk::Math::abs((4 + noProj) % 8 - 4.) - 2.); - center[1] = 0.; - center[2] = 0.; - e2->SetInput(e1->GetOutput()); - e2->SetGeometry(oneProjGeometry); - e2->SetDensity(-1.); - e2->SetAxis(semiprincipalaxis); - e2->SetCenter(center); - e2->SetAngle(0.); - e2->Update(); - - // Adding each projection to the projection stack - if (noProj > 0) // After the first projection, we use the output as input - { - ProjectionStackType::Pointer wholeImage = pasteFilter->GetOutput(); - wholeImage->DisconnectPipeline(); - pasteFilter->SetDestinationImage(wholeImage); - } - pasteFilter->SetSourceImage(e2->GetOutput()); - pasteFilter->SetSourceRegion(e2->GetOutput()->GetLargestPossibleRegion()); - pasteFilter->SetDestinationIndex(destinationIndex); - pasteFilter->UpdateLargestPossibleRegion(); - destinationIndex[2]++; - - // Signal - signalFile << (noProj % 8) / 8. << std::endl; - } - signalFile.close(); - - // Create a vector field and its (very rough) inverse - using IteratorType = itk::ImageRegionIteratorWithIndex; - - auto deformationField = DVFSequenceImageType::New(); - auto inverseDeformationField = DVFSequenceImageType::New(); - - auto sizeMotion = itk::MakeSize(fourDSize[0], fourDSize[1], fourDSize[2], 2); - auto originMotion = itk::MakePoint(-63., -31., -63., 0.); - DVFSequenceImageType::RegionType regionMotion; - regionMotion.SetSize(sizeMotion); - - auto spacingMotion = itk::MakeVector(fourDSpacing[0], fourDSpacing[1], fourDSpacing[2], fourDSpacing[3]); - deformationField->SetRegions(regionMotion); - deformationField->SetOrigin(originMotion); - deformationField->SetSpacing(spacingMotion); - deformationField->Allocate(); - - inverseDeformationField->SetRegions(regionMotion); - inverseDeformationField->SetOrigin(originMotion); - inverseDeformationField->SetSpacing(spacingMotion); - inverseDeformationField->Allocate(); - - // Vector Field initilization - DVFVectorType vec; - IteratorType dvfIt(deformationField, deformationField->GetLargestPossibleRegion()); - IteratorType idvfIt(inverseDeformationField, inverseDeformationField->GetLargestPossibleRegion()); - - DVFSequenceImageType::OffsetType DVFCenter; - DVFSequenceImageType::IndexType toCenter; - DVFCenter.Fill(0); - DVFCenter[0] = sizeMotion[0] / 2; - DVFCenter[1] = sizeMotion[1] / 2; - DVFCenter[2] = sizeMotion[2] / 2; - while (!dvfIt.IsAtEnd()) - { - vec.Fill(0.); - toCenter = dvfIt.GetIndex() - DVFCenter; - - if (0.3 * toCenter[0] * toCenter[0] + toCenter[1] * toCenter[1] + toCenter[2] * toCenter[2] < 40) - { - if (dvfIt.GetIndex()[3] == 0) - vec[0] = -8.; - else - vec[0] = 8.; - } - dvfIt.Set(vec); - idvfIt.Set(-vec); - - ++dvfIt; - ++idvfIt; - } - - // Ground truth - auto * Volumes = new VolumeType::Pointer[fourDSize[3]]; - auto join = itk::JoinSeriesImageFilter::New(); - - for (itk::SizeValueType n = 0; n < fourDSize[3]; n++) - { - // Ellipse 1 - using DEType = rtk::DrawEllipsoidImageFilter; - auto de1 = DEType::New(); - de1->SetInput(tomographySource->GetOutput()); - de1->SetDensity(2.); - DEType::VectorType axis; - axis.Fill(60.); - axis[1] = 30; - de1->SetAxis(axis); - DEType::VectorType center; - center.Fill(0.); - de1->SetCenter(center); - de1->SetAngle(0.); - de1->InPlaceOff(); - TRY_AND_EXIT_ON_ITK_EXCEPTION(de1->Update()) - - // Ellipse 2 - auto de2 = DEType::New(); - de2->SetInput(de1->GetOutput()); - de2->SetDensity(-1.); - DEType::VectorType axis2; - axis2.Fill(8.); - de2->SetAxis(axis2); - DEType::VectorType center2; - center2[0] = 4 * (itk::Math::abs((4 + n) % 8 - 4.) - 2.); - center2[1] = 0.; - center2[2] = 0.; - de2->SetCenter(center2); - de2->SetAngle(0.); - de2->InPlaceOff(); - TRY_AND_EXIT_ON_ITK_EXCEPTION(de2->Update()); - - Volumes[n] = de2->GetOutput(); - Volumes[n]->DisconnectPipeline(); - join->SetInput(n, Volumes[n]); - } - join->Update(); + auto data = rtk::GenerateFourDTestData(FAST_TESTS_NO_CHECKS); // ROI using DEType = rtk::DrawEllipsoidImageFilter; auto roi = DEType::New(); - roi->SetInput(tomographySource->GetOutput()); + roi->SetInput(data.SingleVolume); roi->SetDensity(1.); DEType::VectorType axis; axis.Fill(15.); @@ -285,19 +57,18 @@ main(int, char **) // Read the phases file auto phaseReader = rtk::PhasesToInterpolationWeights::New(); - phaseReader->SetFileName(signalFileName); - phaseReader->SetNumberOfReconstructedFrames(fourDSize[3]); + phaseReader->SetFileName(data.SignalFileName); + phaseReader->SetNumberOfReconstructedFrames(data.InitialVolumeSeries->GetLargestPossibleRegion().GetSize(3)); phaseReader->Update(); // Set the forward and back projection filters to be used using ROOSTERFilterType = rtk::FourDROOSTERConeBeamReconstructionFilter; auto rooster = ROOSTERFilterType::New(); - rooster->SetInputVolumeSeries(fourdSource->GetOutput()); - rooster->SetInputProjectionStack(pasteFilter->GetOutput()); - rooster->SetGeometry(geometry); + rooster->SetInputVolumeSeries(data.InitialVolumeSeries); + rooster->SetInputProjectionStack(data.Projections); + rooster->SetGeometry(data.Geometry); rooster->SetWeights(phaseReader->GetOutput()); - rooster->SetSignal(rtk::ReadSignalFile(signalFileName)); - rooster->SetGeometry(geometry); + rooster->SetSignal(rtk::ReadSignalFile(data.SignalFileName)); rooster->SetCG_iterations(2); rooster->SetMainLoop_iterations(2); @@ -331,7 +102,7 @@ main(int, char **) TRY_AND_EXIT_ON_ITK_EXCEPTION(rooster->Update()); - CheckImageQuality(rooster->GetOutput(), join->GetOutput(), 0.25, 15, 2.0); + CheckImageQuality(rooster->GetOutput(), data.GroundTruth, 0.25, 15, 2.0); std::cout << "\n\nTest PASSED! " << std::endl; std::cout << "\n\n****** Case 2: Joseph forward projector, voxel-based back projector, positivity, no motion mask, " @@ -349,13 +120,13 @@ main(int, char **) rooster->SetPerformTVTemporalDenoising(false); rooster->SetPerformL0TemporalDenoising(true); rooster->SetPerformWarping(true); - rooster->SetDisplacementField(deformationField); + rooster->SetDisplacementField(data.DVF); rooster->SetComputeInverseWarpingByConjugateGradient(true); rooster->SetUseNearestNeighborInterpolationInWarping(true); TRY_AND_EXIT_ON_ITK_EXCEPTION(rooster->Update()); - CheckImageQuality(rooster->GetOutput(), join->GetOutput(), 0.25, 15, 2.0); + CheckImageQuality(rooster->GetOutput(), data.GroundTruth, 0.25, 15, 2.0); std::cout << "\n\nTest PASSED! " << std::endl; std::cout << "\n\n****** Case 3: Joseph forward projector, voxel-based back projector, no positivity, motion mask, " @@ -374,14 +145,14 @@ main(int, char **) rooster->SetPerformTVTemporalDenoising(true); rooster->SetPerformL0TemporalDenoising(false); rooster->SetPerformWarping(true); - rooster->SetDisplacementField(deformationField); + rooster->SetDisplacementField(data.DVF); rooster->SetComputeInverseWarpingByConjugateGradient(false); - rooster->SetInverseDisplacementField(inverseDeformationField); + rooster->SetInverseDisplacementField(data.InverseDVF); rooster->SetUseNearestNeighborInterpolationInWarping(false); TRY_AND_EXIT_ON_ITK_EXCEPTION(rooster->Update()); - CheckImageQuality(rooster->GetOutput(), join->GetOutput(), 0.25, 15, 2.0); + CheckImageQuality(rooster->GetOutput(), data.GroundTruth, 0.25, 15, 2.0); std::cout << "\n\nTest PASSED! " << std::endl; #ifdef USE_CUDA @@ -402,12 +173,11 @@ main(int, char **) TRY_AND_EXIT_ON_ITK_EXCEPTION(rooster->Update()); - CheckImageQuality(rooster->GetOutput(), join->GetOutput(), 0.25, 15, 2.0); + CheckImageQuality(rooster->GetOutput(), data.GroundTruth, 0.25, 15, 2.0); std::cout << "\n\nTest PASSED! " << std::endl; #endif - itksys::SystemTools::RemoveFile(signalFileName.c_str()); - delete[] Volumes; + itksys::SystemTools::RemoveFile(data.SignalFileName.c_str()); return EXIT_SUCCESS; } From ead2406dd7e217000ca9402022267d15616f5528 Mon Sep 17 00:00:00 2001 From: Cyril Mory Date: Thu, 8 Jan 2026 13:44:40 +0100 Subject: [PATCH 2/3] ENH: Add FourDConjugateGradient example and corresponding doc page --- .../rtkfourdconjugategradient.cxx | 1 - .../FourDConjugateGradient/CMakeLists.txt | 12 +++ .../FourDConjugateGradient.cxx | 79 +++++++++++++++++++ examples/FourDConjugateGradient/README.md | 17 ++++ test/CMakeLists.txt | 3 + 5 files changed, 111 insertions(+), 1 deletion(-) create mode 100644 examples/FourDConjugateGradient/CMakeLists.txt create mode 100644 examples/FourDConjugateGradient/FourDConjugateGradient.cxx create mode 100644 examples/FourDConjugateGradient/README.md diff --git a/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx b/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx index d181ee40d..aeedfbe3a 100644 --- a/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx +++ b/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx @@ -26,7 +26,6 @@ #ifdef RTK_USE_CUDA # include "itkCudaImage.h" -# include "rtkCudaConstantVolumeSeriesSource.h" #endif #include diff --git a/examples/FourDConjugateGradient/CMakeLists.txt b/examples/FourDConjugateGradient/CMakeLists.txt new file mode 100644 index 000000000..1941d871d --- /dev/null +++ b/examples/FourDConjugateGradient/CMakeLists.txt @@ -0,0 +1,12 @@ +cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR) + +# This project is designed to be built outside the RTK source tree. +project(FourDConjugateGradient) + +# Find ITK with RTK +find_package(ITK REQUIRED COMPONENTS RTK) +include(${ITK_USE_FILE}) + +# Executable(s) +add_executable(FourDConjugateGradient FourDConjugateGradient.cxx) +target_link_libraries(FourDConjugateGradient ${ITK_LIBRARIES}) diff --git a/examples/FourDConjugateGradient/FourDConjugateGradient.cxx b/examples/FourDConjugateGradient/FourDConjugateGradient.cxx new file mode 100644 index 000000000..b34854780 --- /dev/null +++ b/examples/FourDConjugateGradient/FourDConjugateGradient.cxx @@ -0,0 +1,79 @@ +// #include "rtkThreeDCircularProjectionGeometryXMLFile.h" +#include "rtkFourDConjugateGradientConeBeamReconstructionFilter.h" +#include "rtkIterationCommands.h" +#include "rtkSignalToInterpolationWeights.h" +#include "rtkReorderProjectionsImageFilter.h" +#include "../../test/rtkFourDTestHelper.h" + +#ifdef RTK_USE_CUDA +# include +#endif +#include + +int +main(int argc, char * argv[]) +{ + using OutputPixelType = float; + +#ifdef RTK_USE_CUDA + using VolumeSeriesType = itk::CudaImage; + using ProjectionStackType = itk::CudaImage; + using VolumeType = itk::CudaImage; +#else + using VolumeSeriesType = itk::Image; + using ProjectionStackType = itk::Image; + using VolumeType = itk::Image; +#endif + + auto data = rtk::GenerateFourDTestData(FAST_TESTS_NO_CHECKS); + + // Re-order geometry and projections + // In the new order, projections with identical phases are packed together + auto reorder = rtk::ReorderProjectionsImageFilter::New(); + reorder->SetInput(data.Projections); + reorder->SetInputGeometry(data.Geometry); + reorder->SetInputSignal(data.Signal); + TRY_AND_EXIT_ON_ITK_EXCEPTION(reorder->Update()) + + // Release the memory holding the stack of original projections + data.Projections->ReleaseData(); + + // Compute the interpolation weights + auto signalToInterpolationWeights = rtk::SignalToInterpolationWeights::New(); + signalToInterpolationWeights->SetSignal(reorder->GetOutputSignal()); + signalToInterpolationWeights->SetNumberOfReconstructedFrames( + data.InitialVolumeSeries->GetLargestPossibleRegion().GetSize(3)); + TRY_AND_EXIT_ON_ITK_EXCEPTION(signalToInterpolationWeights->Update()) + + // Set the forward and back projection filters to be used + using ConjugateGradientFilterType = + rtk::FourDConjugateGradientConeBeamReconstructionFilter; + auto fourdconjugategradient = ConjugateGradientFilterType::New(); + + fourdconjugategradient->SetInputVolumeSeries(data.InitialVolumeSeries); + fourdconjugategradient->SetNumberOfIterations(2); +#ifdef RTK_USE_CUDA + fourdconjugategradient->SetCudaConjugateGradient(true); + fourdconjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_CUDARAYCAST); + fourdconjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_CUDAVOXELBASED); +#else + fourdconjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_JOSEPH); + fourdconjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_VOXELBASED); +#endif + + // Set the newly ordered arguments + fourdconjugategradient->SetInputProjectionStack(reorder->GetOutput()); + fourdconjugategradient->SetGeometry(reorder->GetOutputGeometry()); + fourdconjugategradient->SetWeights(signalToInterpolationWeights->GetOutput()); + fourdconjugategradient->SetSignal(reorder->GetOutputSignal()); + + auto verboseIterationCommand = rtk::VerboseIterationCommand::New(); + fourdconjugategradient->AddObserver(itk::AnyEvent(), verboseIterationCommand); + + TRY_AND_EXIT_ON_ITK_EXCEPTION(fourdconjugategradient->Update()) + + // Write + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(fourdconjugategradient->GetOutput(), "fourdconjugategradient.mha")); + + return EXIT_SUCCESS; +} diff --git a/examples/FourDConjugateGradient/README.md b/examples/FourDConjugateGradient/README.md new file mode 100644 index 000000000..2df9c10f9 --- /dev/null +++ b/examples/FourDConjugateGradient/README.md @@ -0,0 +1,17 @@ +# 4D Conjugate Gradient + +FourDConjugateGradient shows how to perform iterative cone-beam CT reconstruction using either CPU or GPU resources. +You can refer to the [projectors documentation](../../documentation/docs/Projectors.md) to see all options available for the back and forwardprojections. + +This example generates its own input data: a phantom made of two ellipsoids, one of which is moving. Projections are computed through this moving phantom. + +`````{tab-set} + +````{tab-item} C++ + +```{literalinclude} ./FourDConjugateGradient.cxx +:language: c++ +``` +```` + +````` diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 107d7016a..67a084c10 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -376,6 +376,9 @@ rtk_add_test(rtkConjugateGradientExampleTest ../examples/ConjugateGradient/ConjugateGradient.cxx DATA{Input/Forbild/Thorax} ) +rtk_add_test(rtkFourDConjugateGradientExampleTest + ../examples/FourDConjugateGradient/FourDConjugateGradient.cxx +) if(ITK_WRAP_PYTHON) itk_python_add_test(NAME rtkMaximumIntensityPythonTest COMMAND rtkMaximumIntensity.py) From b928a5f17cd9b416ea4ce165afdbb0260cdce7fc Mon Sep 17 00:00:00 2001 From: Cyril Mory Date: Mon, 19 Jan 2026 07:26:22 +0100 Subject: [PATCH 3/3] ENH: Move 4D example data generation to a separate example Adapt FourDConjugateGradient example to use the example 4D data Order the tests so that GenerateFourDData runs first. It saves its data in the test/ folder, where tests read --- .../FourDConjugateGradient.cxx | 45 +++- examples/GenerateFourDData/CMakeLists.txt | 12 + .../GenerateFourDData/GenerateFourDData.cxx | 228 ++++++++++++++++++ examples/GenerateFourDData/README.md | 21 ++ test/CMakeLists.txt | 15 ++ 5 files changed, 310 insertions(+), 11 deletions(-) create mode 100644 examples/GenerateFourDData/CMakeLists.txt create mode 100644 examples/GenerateFourDData/GenerateFourDData.cxx create mode 100644 examples/GenerateFourDData/README.md diff --git a/examples/FourDConjugateGradient/FourDConjugateGradient.cxx b/examples/FourDConjugateGradient/FourDConjugateGradient.cxx index b34854780..e36e04f1b 100644 --- a/examples/FourDConjugateGradient/FourDConjugateGradient.cxx +++ b/examples/FourDConjugateGradient/FourDConjugateGradient.cxx @@ -1,9 +1,9 @@ -// #include "rtkThreeDCircularProjectionGeometryXMLFile.h" #include "rtkFourDConjugateGradientConeBeamReconstructionFilter.h" #include "rtkIterationCommands.h" #include "rtkSignalToInterpolationWeights.h" #include "rtkReorderProjectionsImageFilter.h" -#include "../../test/rtkFourDTestHelper.h" +#include "rtkProjectionsReader.h" +#include "rtkThreeDCircularProjectionGeometryXMLFileReader.h" #ifdef RTK_USE_CUDA # include @@ -25,24 +25,47 @@ main(int argc, char * argv[]) using VolumeType = itk::Image; #endif - auto data = rtk::GenerateFourDTestData(FAST_TESTS_NO_CHECKS); + // Generate the input volume series, used as initial estimate by 4D conjugate gradient + auto fourDOrigin = itk::MakePoint(-63., -31., -63., 0.); + auto fourDSpacing = itk::MakeVector(4., 4., 4., 1.); + auto fourDSize = itk::MakeSize(32, 16, 32, 8); + + using ConstantFourDSourceType = rtk::ConstantImageSource; + auto fourDSource = ConstantFourDSourceType::New(); + + fourDSource->SetOrigin(fourDOrigin); + fourDSource->SetSpacing(fourDSpacing); + fourDSource->SetSize(fourDSize); + fourDSource->SetConstant(0.); + fourDSource->Update(); + + // Read geometry, projections and signal + rtk::ThreeDCircularProjectionGeometry::Pointer geometry; + TRY_AND_EXIT_ON_ITK_EXCEPTION(geometry = rtk::ReadGeometry("four_d_geometry.xml")); + + using ReaderType = rtk::ProjectionsReader; + auto projectionsReader = ReaderType::New(); + std::vector fileNames = std::vector(); + fileNames.push_back("four_d_projections.mha"); + projectionsReader->SetFileNames(fileNames); + TRY_AND_EXIT_ON_ITK_EXCEPTION(projectionsReader->Update()); // Re-order geometry and projections // In the new order, projections with identical phases are packed together - auto reorder = rtk::ReorderProjectionsImageFilter::New(); - reorder->SetInput(data.Projections); - reorder->SetInputGeometry(data.Geometry); - reorder->SetInputSignal(data.Signal); + std::vector signal = rtk::ReadSignalFile("four_d_signal.txt"); + auto reorder = rtk::ReorderProjectionsImageFilter::New(); + reorder->SetInput(projectionsReader->GetOutput()); + reorder->SetInputGeometry(geometry); + reorder->SetInputSignal(signal); TRY_AND_EXIT_ON_ITK_EXCEPTION(reorder->Update()) // Release the memory holding the stack of original projections - data.Projections->ReleaseData(); + projectionsReader->GetOutput()->ReleaseData(); // Compute the interpolation weights auto signalToInterpolationWeights = rtk::SignalToInterpolationWeights::New(); signalToInterpolationWeights->SetSignal(reorder->GetOutputSignal()); - signalToInterpolationWeights->SetNumberOfReconstructedFrames( - data.InitialVolumeSeries->GetLargestPossibleRegion().GetSize(3)); + signalToInterpolationWeights->SetNumberOfReconstructedFrames(fourDSize[3]); TRY_AND_EXIT_ON_ITK_EXCEPTION(signalToInterpolationWeights->Update()) // Set the forward and back projection filters to be used @@ -50,7 +73,7 @@ main(int argc, char * argv[]) rtk::FourDConjugateGradientConeBeamReconstructionFilter; auto fourdconjugategradient = ConjugateGradientFilterType::New(); - fourdconjugategradient->SetInputVolumeSeries(data.InitialVolumeSeries); + fourdconjugategradient->SetInputVolumeSeries(fourDSource->GetOutput()); fourdconjugategradient->SetNumberOfIterations(2); #ifdef RTK_USE_CUDA fourdconjugategradient->SetCudaConjugateGradient(true); diff --git a/examples/GenerateFourDData/CMakeLists.txt b/examples/GenerateFourDData/CMakeLists.txt new file mode 100644 index 000000000..04b65011a --- /dev/null +++ b/examples/GenerateFourDData/CMakeLists.txt @@ -0,0 +1,12 @@ +cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR) + +# This project is designed to be built outside the RTK source tree. +project(GenerateFourDData) + +# Find ITK with RTK +find_package(ITK REQUIRED COMPONENTS RTK) +include(${ITK_USE_FILE}) + +# Executable(s) +add_executable(GenerateFourDData GenerateFourDData.cxx) +target_link_libraries(GenerateFourDData ${ITK_LIBRARIES}) diff --git a/examples/GenerateFourDData/GenerateFourDData.cxx b/examples/GenerateFourDData/GenerateFourDData.cxx new file mode 100644 index 000000000..044f74bcd --- /dev/null +++ b/examples/GenerateFourDData/GenerateFourDData.cxx @@ -0,0 +1,228 @@ +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "rtkConstantImageSource.h" +#include "rtkRayEllipsoidIntersectionImageFilter.h" +#include "rtkDrawEllipsoidImageFilter.h" +#include "rtkThreeDCircularProjectionGeometry.h" +#ifdef RTK_USE_CUDA +# include +#endif + +int +main(int argc, char * argv[]) +{ + using PixelType = float; + using DVFVectorType = itk::CovariantVector; + +#ifdef RTK_USE_CUDA + using VolumeSeriesType = itk::CudaImage; + using ProjectionStackType = itk::CudaImage; + using VolumeType = itk::CudaImage; + using DVFSequenceImageType = itk::CudaImage; +#else + using VolumeSeriesType = itk::Image; + using ProjectionStackType = itk::Image; + using VolumeType = itk::Image; + using DVFSequenceImageType = itk::Image; +#endif + using GeometryType = rtk::ThreeDCircularProjectionGeometry; + constexpr unsigned int NumberOfProjectionImages = 64; + + /* ========================= + * Constant volume source + * ========================= */ + using ConstantVolumeSourceType = rtk::ConstantImageSource; + auto volumeSource = ConstantVolumeSourceType::New(); + + auto origin = itk::MakePoint(-63., -31., -63.); + auto spacing = itk::MakeVector(4., 4., 4.); + auto size = itk::MakeSize(32, 16, 32); + + volumeSource->SetOrigin(origin); + volumeSource->SetSpacing(spacing); + volumeSource->SetSize(size); + volumeSource->SetConstant(0.); + volumeSource->Update(); + + /* ========================= + * Projection accumulation + * ========================= */ + using ConstantProjectionSourceType = rtk::ConstantImageSource; + auto projectionsSource = ConstantProjectionSourceType::New(); + + auto projOrigin = itk::MakePoint(-254., -254., -254.); + auto projSpacing = itk::MakeVector(8., 8., 1.); + auto projSize = itk::MakeSize(64, 64, NumberOfProjectionImages); + + projectionsSource->SetOrigin(projOrigin); + projectionsSource->SetSpacing(projSpacing); + projectionsSource->SetSize(projSize); + projectionsSource->SetConstant(0.); + projectionsSource->Update(); + + auto oneProjectionSource = ConstantProjectionSourceType::New(); + projSize[2] = 1; + oneProjectionSource->SetOrigin(projOrigin); + oneProjectionSource->SetSpacing(projSpacing); + oneProjectionSource->SetSize(projSize); + oneProjectionSource->SetConstant(0.); + + using REIType = rtk::RayEllipsoidIntersectionImageFilter; + using PasteType = itk::PasteImageFilter; + + /* Allocate explicit accumulator image */ + auto accumulated = ProjectionStackType::New(); + accumulated->SetOrigin(projectionsSource->GetOutput()->GetOrigin()); + accumulated->SetSpacing(projectionsSource->GetOutput()->GetSpacing()); + accumulated->SetDirection(projectionsSource->GetOutput()->GetDirection()); + accumulated->SetRegions(projectionsSource->GetOutput()->GetLargestPossibleRegion()); + accumulated->Allocate(); + accumulated->FillBuffer(0.); + + auto paste = PasteType::New(); + paste->SetDestinationImage(accumulated); + + itk::Index<3> destIndex; + destIndex.Fill(0); + + // Create the signal file and the geometry, to be filled in the for loop + std::string signalFileName = "four_d_signal.txt"; + std::ofstream signalFile(signalFileName.c_str()); + auto geometry = GeometryType::New(); + + for (unsigned int i = 0; i < NumberOfProjectionImages; ++i) + { + geometry->AddProjection(600., 1200., i * 360. / NumberOfProjectionImages, 0, 0, 0, 0, 20, 15); + + auto geom = GeometryType::New(); + geom->AddProjection(600., 1200., i * 360. / NumberOfProjectionImages, 0, 0, 0, 0, 20, 15); + + auto e1 = REIType::New(); + e1->SetInput(oneProjectionSource->GetOutput()); + e1->SetGeometry(geom); + e1->SetDensity(2.); + e1->SetAxis(itk::MakeVector(60., 30., 60.)); + e1->SetCenter(itk::MakeVector(0., 0., 0.)); + e1->InPlaceOff(); + e1->Update(); + + auto e2 = REIType::New(); + e2->SetInput(e1->GetOutput()); + e2->SetGeometry(geom); + e2->SetDensity(-1.); + e2->SetAxis(itk::MakeVector(8., 8., 8.)); + auto center = itk::MakeVector(4 * (itk::Math::abs((4 + i) % 8 - 4.) - 2.), 0., 0.); + e2->SetCenter(center); + e2->InPlaceOff(); + e2->Update(); + + paste->SetSourceImage(e2->GetOutput()); + paste->SetSourceRegion(e2->GetOutput()->GetLargestPossibleRegion()); + paste->SetDestinationIndex(destIndex); + paste->Update(); + + accumulated = paste->GetOutput(); + paste->SetDestinationImage(accumulated); + + destIndex[2]++; + + signalFile << (i % 8) / 8. << std::endl; + } + + signalFile.close(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(accumulated, "four_d_projections.mha")); + TRY_AND_EXIT_ON_ITK_EXCEPTION(rtk::WriteGeometry(geometry, "four_d_geometry.xml")) + + /* ========================= + * DVF & inverse DVF + * ========================= */ + auto fourDOrigin = itk::MakePoint(-63., -31., -63., 0.); + auto fourDSpacing = itk::MakeVector(4., 4., 4., 1.); + auto fourDSize = itk::MakeSize(32, 16, 32, 8); + + auto dvf = DVFSequenceImageType::New(); + auto idvf = DVFSequenceImageType::New(); + + typename DVFSequenceImageType::RegionType region; + auto dvfSize = itk::MakeSize(fourDSize[0], fourDSize[1], fourDSize[2], 2); + region.SetSize(dvfSize); + + dvf->SetRegions(region); + dvf->SetOrigin(fourDOrigin); + dvf->SetSpacing(fourDSpacing); + dvf->Allocate(); + + idvf->SetRegions(region); + idvf->SetOrigin(fourDOrigin); + idvf->SetSpacing(fourDSpacing); + idvf->Allocate(); + + itk::ImageRegionIteratorWithIndex it(dvf, region); + itk::ImageRegionIteratorWithIndex iit(idvf, region); + + DVFVectorType v; + typename DVFSequenceImageType::IndexType centerIndex; + centerIndex.Fill(0); + centerIndex[0] = dvfSize[0] / 2; + centerIndex[1] = dvfSize[1] / 2; + centerIndex[2] = dvfSize[2] / 2; + + for (; !it.IsAtEnd(); ++it, ++iit) + { + v.Fill(0.); + auto d = it.GetIndex() - centerIndex; + if (0.3 * d[0] * d[0] + d[1] * d[1] + d[2] * d[2] < 40) + v[0] = (it.GetIndex()[3] == 0 ? -8. : 8.); + it.Set(v); + iit.Set(-v); + } + + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(dvf, "four_d_dvf.mha")); + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(idvf, "four_d_idvf.mha")); + + /* ========================= + * Ground truth + * ========================= */ + auto join = itk::JoinSeriesImageFilter::New(); + + for (unsigned int t = 0; t < fourDSize[3]; ++t) + { + using DEType = rtk::DrawEllipsoidImageFilter; + + auto de1 = DEType::New(); + de1->SetInput(volumeSource->GetOutput()); + de1->SetDensity(2.); + de1->SetAxis(itk::MakeVector(60., 30., 60.)); + de1->SetCenter(itk::MakeVector(0., 0., 0.)); + de1->InPlaceOff(); + de1->Update(); + + auto de2 = DEType::New(); + de2->SetInput(de1->GetOutput()); + de2->SetDensity(-1.); + de2->SetAxis(itk::MakeVector(8., 8., 8.)); + de2->SetCenter(itk::MakeVector(4 * (itk::Math::abs((4 + t) % 8 - 4.) - 2.), 0., 0.)); + de2->InPlaceOff(); + de2->Update(); + + using DuplicatorType = itk::ImageDuplicator; + auto duplicator = DuplicatorType::New(); + duplicator->SetInputImage(de2->GetOutput()); + duplicator->Update(); + + join->SetInput(t, duplicator->GetOutput()); + } + + join->Update(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(join->GetOutput(), "four_d_ground_truth.mha")); + + return EXIT_SUCCESS; +} diff --git a/examples/GenerateFourDData/README.md b/examples/GenerateFourDData/README.md new file mode 100644 index 000000000..c416c1056 --- /dev/null +++ b/examples/GenerateFourDData/README.md @@ -0,0 +1,21 @@ +# Generate 4D data + +GenerateFourDData shows how to generate a full set of input data to be used in the examples of 4D reconstruction methods. It contains: +- a moving phantom, made of two ellipsoids, one of which moves along the first dimension +- a geometry +- a phase signal +- the set of projections computed through the moving phantom +- a deformation vector field describing the motion of the moving ellipsoid +- the inverse deformation vector field + +You can also skip this part and [download the data](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531). + +`````{tab-set} + +````{tab-item} C++ + +```{literalinclude} ./GenerateFourDData.cxx +:language: c++ +``` +```` +````` diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 67a084c10..32e66367b 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -376,9 +376,24 @@ rtk_add_test(rtkConjugateGradientExampleTest ../examples/ConjugateGradient/ConjugateGradient.cxx DATA{Input/Forbild/Thorax} ) +rtk_add_test(rtkGenerateFourDDataExampleTest + ../examples/GenerateFourDData/GenerateFourDData.cxx +) +set_tests_properties( + rtkGenerateFourDDataExampleTest + PROPERTIES + FIXTURES_SETUP + FourDData +) rtk_add_test(rtkFourDConjugateGradientExampleTest ../examples/FourDConjugateGradient/FourDConjugateGradient.cxx ) +set_tests_properties( + rtkFourDConjugateGradientExampleTest + PROPERTIES + FIXTURES_REQUIRED + FourDData +) if(ITK_WRAP_PYTHON) itk_python_add_test(NAME rtkMaximumIntensityPythonTest COMMAND rtkMaximumIntensity.py)