Skip to content

Commit b928a5f

Browse files
cyrilmorycyrilmory
committed
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
1 parent ead2406 commit b928a5f

File tree

5 files changed

+310
-11
lines changed

5 files changed

+310
-11
lines changed

‎examples/FourDConjugateGradient/FourDConjugateGradient.cxx‎

Lines changed: 34 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,9 @@
1-
// #include "rtkThreeDCircularProjectionGeometryXMLFile.h"
21
#include "rtkFourDConjugateGradientConeBeamReconstructionFilter.h"
32
#include "rtkIterationCommands.h"
43
#include "rtkSignalToInterpolationWeights.h"
54
#include "rtkReorderProjectionsImageFilter.h"
6-
#include "../../test/rtkFourDTestHelper.h"
5+
#include "rtkProjectionsReader.h"
6+
#include "rtkThreeDCircularProjectionGeometryXMLFileReader.h"
77

88
#ifdef RTK_USE_CUDA
99
# include <itkCudaImage.h>
@@ -25,32 +25,55 @@ main(int argc, char * argv[])
2525
using VolumeType = itk::Image<OutputPixelType, 3>;
2626
#endif
2727

28-
auto data = rtk::GenerateFourDTestData<OutputPixelType>(FAST_TESTS_NO_CHECKS);
28+
// Generate the input volume series, used as initial estimate by 4D conjugate gradient
29+
auto fourDOrigin = itk::MakePoint(-63., -31., -63., 0.);
30+
auto fourDSpacing = itk::MakeVector(4., 4., 4., 1.);
31+
auto fourDSize = itk::MakeSize(32, 16, 32, 8);
32+
33+
using ConstantFourDSourceType = rtk::ConstantImageSource<VolumeSeriesType>;
34+
auto fourDSource = ConstantFourDSourceType::New();
35+
36+
fourDSource->SetOrigin(fourDOrigin);
37+
fourDSource->SetSpacing(fourDSpacing);
38+
fourDSource->SetSize(fourDSize);
39+
fourDSource->SetConstant(0.);
40+
fourDSource->Update();
41+
42+
// Read geometry, projections and signal
43+
rtk::ThreeDCircularProjectionGeometry::Pointer geometry;
44+
TRY_AND_EXIT_ON_ITK_EXCEPTION(geometry = rtk::ReadGeometry("four_d_geometry.xml"));
45+
46+
using ReaderType = rtk::ProjectionsReader<ProjectionStackType>;
47+
auto projectionsReader = ReaderType::New();
48+
std::vector<std::string> fileNames = std::vector<std::string>();
49+
fileNames.push_back("four_d_projections.mha");
50+
projectionsReader->SetFileNames(fileNames);
51+
TRY_AND_EXIT_ON_ITK_EXCEPTION(projectionsReader->Update());
2952

3053
// Re-order geometry and projections
3154
// In the new order, projections with identical phases are packed together
32-
auto reorder = rtk::ReorderProjectionsImageFilter<ProjectionStackType>::New();
33-
reorder->SetInput(data.Projections);
34-
reorder->SetInputGeometry(data.Geometry);
35-
reorder->SetInputSignal(data.Signal);
55+
std::vector<double> signal = rtk::ReadSignalFile("four_d_signal.txt");
56+
auto reorder = rtk::ReorderProjectionsImageFilter<ProjectionStackType>::New();
57+
reorder->SetInput(projectionsReader->GetOutput());
58+
reorder->SetInputGeometry(geometry);
59+
reorder->SetInputSignal(signal);
3660
TRY_AND_EXIT_ON_ITK_EXCEPTION(reorder->Update())
3761

3862
// Release the memory holding the stack of original projections
39-
data.Projections->ReleaseData();
63+
projectionsReader->GetOutput()->ReleaseData();
4064

4165
// Compute the interpolation weights
4266
auto signalToInterpolationWeights = rtk::SignalToInterpolationWeights::New();
4367
signalToInterpolationWeights->SetSignal(reorder->GetOutputSignal());
44-
signalToInterpolationWeights->SetNumberOfReconstructedFrames(
45-
data.InitialVolumeSeries->GetLargestPossibleRegion().GetSize(3));
68+
signalToInterpolationWeights->SetNumberOfReconstructedFrames(fourDSize[3]);
4669
TRY_AND_EXIT_ON_ITK_EXCEPTION(signalToInterpolationWeights->Update())
4770

4871
// Set the forward and back projection filters to be used
4972
using ConjugateGradientFilterType =
5073
rtk::FourDConjugateGradientConeBeamReconstructionFilter<VolumeSeriesType, ProjectionStackType>;
5174
auto fourdconjugategradient = ConjugateGradientFilterType::New();
5275

53-
fourdconjugategradient->SetInputVolumeSeries(data.InitialVolumeSeries);
76+
fourdconjugategradient->SetInputVolumeSeries(fourDSource->GetOutput());
5477
fourdconjugategradient->SetNumberOfIterations(2);
5578
#ifdef RTK_USE_CUDA
5679
fourdconjugategradient->SetCudaConjugateGradient(true);

‎examples/GenerateFourDData/CMakeLists.txt‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR)
2+
3+
# This project is designed to be built outside the RTK source tree.
4+
project(GenerateFourDData)
5+
6+
# Find ITK with RTK
7+
find_package(ITK REQUIRED COMPONENTS RTK)
8+
include(${ITK_USE_FILE})
9+
10+
# Executable(s)
11+
add_executable(GenerateFourDData GenerateFourDData.cxx)
12+
target_link_libraries(GenerateFourDData ${ITK_LIBRARIES})

‎examples/GenerateFourDData/GenerateFourDData.cxx‎

Lines changed: 228 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,228 @@
1+
#include <string>
2+
#include <fstream>
3+
4+
#include <itkImage.h>
5+
#include <itkCovariantVector.h>
6+
#include <itkPasteImageFilter.h>
7+
#include <itkJoinSeriesImageFilter.h>
8+
#include <itkImageRegionIteratorWithIndex.h>
9+
#include <itkImageDuplicator.h>
10+
11+
#include "rtkConstantImageSource.h"
12+
#include "rtkRayEllipsoidIntersectionImageFilter.h"
13+
#include "rtkDrawEllipsoidImageFilter.h"
14+
#include "rtkThreeDCircularProjectionGeometry.h"
15+
#ifdef RTK_USE_CUDA
16+
# include <itkCudaImage.h>
17+
#endif
18+
19+
int
20+
main(int argc, char * argv[])
21+
{
22+
using PixelType = float;
23+
using DVFVectorType = itk::CovariantVector<PixelType, 3>;
24+
25+
#ifdef RTK_USE_CUDA
26+
using VolumeSeriesType = itk::CudaImage<PixelType, 4>;
27+
using ProjectionStackType = itk::CudaImage<PixelType, 3>;
28+
using VolumeType = itk::CudaImage<PixelType, 3>;
29+
using DVFSequenceImageType = itk::CudaImage<DVFVectorType, 4>;
30+
#else
31+
using VolumeSeriesType = itk::Image<PixelType, 4>;
32+
using ProjectionStackType = itk::Image<PixelType, 3>;
33+
using VolumeType = itk::Image<PixelType, 3>;
34+
using DVFSequenceImageType = itk::Image<DVFVectorType, 4>;
35+
#endif
36+
using GeometryType = rtk::ThreeDCircularProjectionGeometry;
37+
constexpr unsigned int NumberOfProjectionImages = 64;
38+
39+
/* =========================
40+
* Constant volume source
41+
* ========================= */
42+
using ConstantVolumeSourceType = rtk::ConstantImageSource<VolumeType>;
43+
auto volumeSource = ConstantVolumeSourceType::New();
44+
45+
auto origin = itk::MakePoint(-63., -31., -63.);
46+
auto spacing = itk::MakeVector(4., 4., 4.);
47+
auto size = itk::MakeSize(32, 16, 32);
48+
49+
volumeSource->SetOrigin(origin);
50+
volumeSource->SetSpacing(spacing);
51+
volumeSource->SetSize(size);
52+
volumeSource->SetConstant(0.);
53+
volumeSource->Update();
54+
55+
/* =========================
56+
* Projection accumulation
57+
* ========================= */
58+
using ConstantProjectionSourceType = rtk::ConstantImageSource<ProjectionStackType>;
59+
auto projectionsSource = ConstantProjectionSourceType::New();
60+
61+
auto projOrigin = itk::MakePoint(-254., -254., -254.);
62+
auto projSpacing = itk::MakeVector(8., 8., 1.);
63+
auto projSize = itk::MakeSize(64, 64, NumberOfProjectionImages);
64+
65+
projectionsSource->SetOrigin(projOrigin);
66+
projectionsSource->SetSpacing(projSpacing);
67+
projectionsSource->SetSize(projSize);
68+
projectionsSource->SetConstant(0.);
69+
projectionsSource->Update();
70+
71+
auto oneProjectionSource = ConstantProjectionSourceType::New();
72+
projSize[2] = 1;
73+
oneProjectionSource->SetOrigin(projOrigin);
74+
oneProjectionSource->SetSpacing(projSpacing);
75+
oneProjectionSource->SetSize(projSize);
76+
oneProjectionSource->SetConstant(0.);
77+
78+
using REIType = rtk::RayEllipsoidIntersectionImageFilter<VolumeType, ProjectionStackType>;
79+
using PasteType = itk::PasteImageFilter<ProjectionStackType, ProjectionStackType, ProjectionStackType>;
80+
81+
/* Allocate explicit accumulator image */
82+
auto accumulated = ProjectionStackType::New();
83+
accumulated->SetOrigin(projectionsSource->GetOutput()->GetOrigin());
84+
accumulated->SetSpacing(projectionsSource->GetOutput()->GetSpacing());
85+
accumulated->SetDirection(projectionsSource->GetOutput()->GetDirection());
86+
accumulated->SetRegions(projectionsSource->GetOutput()->GetLargestPossibleRegion());
87+
accumulated->Allocate();
88+
accumulated->FillBuffer(0.);
89+
90+
auto paste = PasteType::New();
91+
paste->SetDestinationImage(accumulated);
92+
93+
itk::Index<3> destIndex;
94+
destIndex.Fill(0);
95+
96+
// Create the signal file and the geometry, to be filled in the for loop
97+
std::string signalFileName = "four_d_signal.txt";
98+
std::ofstream signalFile(signalFileName.c_str());
99+
auto geometry = GeometryType::New();
100+
101+
for (unsigned int i = 0; i < NumberOfProjectionImages; ++i)
102+
{
103+
geometry->AddProjection(600., 1200., i * 360. / NumberOfProjectionImages, 0, 0, 0, 0, 20, 15);
104+
105+
auto geom = GeometryType::New();
106+
geom->AddProjection(600., 1200., i * 360. / NumberOfProjectionImages, 0, 0, 0, 0, 20, 15);
107+
108+
auto e1 = REIType::New();
109+
e1->SetInput(oneProjectionSource->GetOutput());
110+
e1->SetGeometry(geom);
111+
e1->SetDensity(2.);
112+
e1->SetAxis(itk::MakeVector(60., 30., 60.));
113+
e1->SetCenter(itk::MakeVector(0., 0., 0.));
114+
e1->InPlaceOff();
115+
e1->Update();
116+
117+
auto e2 = REIType::New();
118+
e2->SetInput(e1->GetOutput());
119+
e2->SetGeometry(geom);
120+
e2->SetDensity(-1.);
121+
e2->SetAxis(itk::MakeVector(8., 8., 8.));
122+
auto center = itk::MakeVector(4 * (itk::Math::abs((4 + i) % 8 - 4.) - 2.), 0., 0.);
123+
e2->SetCenter(center);
124+
e2->InPlaceOff();
125+
e2->Update();
126+
127+
paste->SetSourceImage(e2->GetOutput());
128+
paste->SetSourceRegion(e2->GetOutput()->GetLargestPossibleRegion());
129+
paste->SetDestinationIndex(destIndex);
130+
paste->Update();
131+
132+
accumulated = paste->GetOutput();
133+
paste->SetDestinationImage(accumulated);
134+
135+
destIndex[2]++;
136+
137+
signalFile << (i % 8) / 8. << std::endl;
138+
}
139+
140+
signalFile.close();
141+
TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(accumulated, "four_d_projections.mha"));
142+
TRY_AND_EXIT_ON_ITK_EXCEPTION(rtk::WriteGeometry(geometry, "four_d_geometry.xml"))
143+
144+
/* =========================
145+
* DVF & inverse DVF
146+
* ========================= */
147+
auto fourDOrigin = itk::MakePoint(-63., -31., -63., 0.);
148+
auto fourDSpacing = itk::MakeVector(4., 4., 4., 1.);
149+
auto fourDSize = itk::MakeSize(32, 16, 32, 8);
150+
151+
auto dvf = DVFSequenceImageType::New();
152+
auto idvf = DVFSequenceImageType::New();
153+
154+
typename DVFSequenceImageType::RegionType region;
155+
auto dvfSize = itk::MakeSize(fourDSize[0], fourDSize[1], fourDSize[2], 2);
156+
region.SetSize(dvfSize);
157+
158+
dvf->SetRegions(region);
159+
dvf->SetOrigin(fourDOrigin);
160+
dvf->SetSpacing(fourDSpacing);
161+
dvf->Allocate();
162+
163+
idvf->SetRegions(region);
164+
idvf->SetOrigin(fourDOrigin);
165+
idvf->SetSpacing(fourDSpacing);
166+
idvf->Allocate();
167+
168+
itk::ImageRegionIteratorWithIndex<DVFSequenceImageType> it(dvf, region);
169+
itk::ImageRegionIteratorWithIndex<DVFSequenceImageType> iit(idvf, region);
170+
171+
DVFVectorType v;
172+
typename DVFSequenceImageType::IndexType centerIndex;
173+
centerIndex.Fill(0);
174+
centerIndex[0] = dvfSize[0] / 2;
175+
centerIndex[1] = dvfSize[1] / 2;
176+
centerIndex[2] = dvfSize[2] / 2;
177+
178+
for (; !it.IsAtEnd(); ++it, ++iit)
179+
{
180+
v.Fill(0.);
181+
auto d = it.GetIndex() - centerIndex;
182+
if (0.3 * d[0] * d[0] + d[1] * d[1] + d[2] * d[2] < 40)
183+
v[0] = (it.GetIndex()[3] == 0 ? -8. : 8.);
184+
it.Set(v);
185+
iit.Set(-v);
186+
}
187+
188+
TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(dvf, "four_d_dvf.mha"));
189+
TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(idvf, "four_d_idvf.mha"));
190+
191+
/* =========================
192+
* Ground truth
193+
* ========================= */
194+
auto join = itk::JoinSeriesImageFilter<VolumeType, VolumeSeriesType>::New();
195+
196+
for (unsigned int t = 0; t < fourDSize[3]; ++t)
197+
{
198+
using DEType = rtk::DrawEllipsoidImageFilter<VolumeType, VolumeType>;
199+
200+
auto de1 = DEType::New();
201+
de1->SetInput(volumeSource->GetOutput());
202+
de1->SetDensity(2.);
203+
de1->SetAxis(itk::MakeVector(60., 30., 60.));
204+
de1->SetCenter(itk::MakeVector(0., 0., 0.));
205+
de1->InPlaceOff();
206+
de1->Update();
207+
208+
auto de2 = DEType::New();
209+
de2->SetInput(de1->GetOutput());
210+
de2->SetDensity(-1.);
211+
de2->SetAxis(itk::MakeVector(8., 8., 8.));
212+
de2->SetCenter(itk::MakeVector(4 * (itk::Math::abs((4 + t) % 8 - 4.) - 2.), 0., 0.));
213+
de2->InPlaceOff();
214+
de2->Update();
215+
216+
using DuplicatorType = itk::ImageDuplicator<VolumeType>;
217+
auto duplicator = DuplicatorType::New();
218+
duplicator->SetInputImage(de2->GetOutput());
219+
duplicator->Update();
220+
221+
join->SetInput(t, duplicator->GetOutput());
222+
}
223+
224+
join->Update();
225+
TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(join->GetOutput(), "four_d_ground_truth.mha"));
226+
227+
return EXIT_SUCCESS;
228+
}

‎examples/GenerateFourDData/README.md‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
# Generate 4D data
2+
3+
GenerateFourDData shows how to generate a full set of input data to be used in the examples of 4D reconstruction methods. It contains:
4+
- a moving phantom, made of two ellipsoids, one of which moves along the first dimension
5+
- a geometry
6+
- a phase signal
7+
- the set of projections computed through the moving phantom
8+
- a deformation vector field describing the motion of the moving ellipsoid
9+
- the inverse deformation vector field
10+
11+
You can also skip this part and [download the data](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531).
12+
13+
`````{tab-set}
14+
15+
````{tab-item} C++
16+
17+
```{literalinclude} ./GenerateFourDData.cxx
18+
:language: c++
19+
```
20+
````
21+
`````

‎test/CMakeLists.txt‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -376,9 +376,24 @@ rtk_add_test(rtkConjugateGradientExampleTest
376376
../examples/ConjugateGradient/ConjugateGradient.cxx
377377
DATA{Input/Forbild/Thorax}
378378
)
379+
rtk_add_test(rtkGenerateFourDDataExampleTest
380+
../examples/GenerateFourDData/GenerateFourDData.cxx
381+
)
382+
set_tests_properties(
383+
rtkGenerateFourDDataExampleTest
384+
PROPERTIES
385+
FIXTURES_SETUP
386+
FourDData
387+
)
379388
rtk_add_test(rtkFourDConjugateGradientExampleTest
380389
../examples/FourDConjugateGradient/FourDConjugateGradient.cxx
381390
)
391+
set_tests_properties(
392+
rtkFourDConjugateGradientExampleTest
393+
PROPERTIES
394+
FIXTURES_REQUIRED
395+
FourDData
396+
)
382397

383398
if(ITK_WRAP_PYTHON)
384399
itk_python_add_test(NAME rtkMaximumIntensityPythonTest COMMAND rtkMaximumIntensity.py)

0 commit comments

Comments
 (0)