Skip to content

Commit 9cd1901

Browse files
thewtexthewtex
authored
Merge pull request #221 from thewtex/resample-label-image-python
DOC: Simplify ResampleSegmentedImage, Python version
2 parents cb1b6e7 + d3ae296 commit 9cd1901

File tree

9 files changed

+169
-176
lines changed

9 files changed

+169
-176
lines changed

‎src/Core/ImageFunction/ResampleSegmentedImage/2th_cthead1.png.sha512‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
c7e0eba27f3c93061503294dd4881111da357dba5e35bba73e0b4114d5f4737ac9493e240ef1dd3b83aa2df33b6828981dcc9878209b79989e3dcbd42469c135

‎src/Core/ImageFunction/ResampleSegmentedImage/CMakeLists.txt‎

Lines changed: 21 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -5,40 +5,36 @@ project(ResampleSegmentedImage)
55
find_package(ITK REQUIRED)
66
include(${ITK_USE_FILE})
77

8-
if(ENABLE_QUICKVIEW)
9-
find_package(VTK REQUIRED
10-
COMPONENTS
11-
vtkRenderingCore
12-
vtkRenderingGL2PSOpenGL2
13-
)
14-
if(VTK_VERSION VERSION_LESS "8.90.0")
15-
include(${VTK_USE_FILE})
16-
endif()
17-
18-
add_executable(ResampleSegmentedImage Code.cxx)
19-
target_link_libraries(ResampleSegmentedImage ${ITK_LIBRARIES} ${VTK_LIBRARIES})
20-
21-
if(NOT VTK_VERSION VERSION_LESS "8.90.0")
22-
vtk_module_autoinit(
23-
TARGETS ResampleSegmentedImage
24-
MODULES ${VTK_LIBRARIES}
25-
)
26-
endif()
27-
else()
28-
add_executable(ResampleSegmentedImage Code.cxx)
29-
target_link_libraries(ResampleSegmentedImage ${ITK_LIBRARIES})
30-
endif()
8+
add_executable(ResampleSegmentedImage Code.cxx)
9+
target_link_libraries(ResampleSegmentedImage ${ITK_LIBRARIES})
3110

3211
install(TARGETS ResampleSegmentedImage
3312
DESTINATION bin/ITKExamples/Core/ImageFunction
3413
COMPONENT Runtime
3514
)
3615

37-
install(FILES Code.cxx CMakeLists.txt
16+
install(FILES Code.cxx CMakeLists.txt Code.py
3817
DESTINATION share/ITKExamples/Code/Core/ImageFunction/ResampleSegmentedImage/
3918
COMPONENT Code
4019
)
4120

4221
enable_testing()
4322
add_test(NAME ResampleSegmentedImageTest
44-
COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/ResampleSegmentedImage Gourds.png 2)
23+
COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/ResampleSegmentedImage
24+
${CMAKE_CURRENT_BINARY_DIR}/2th_cthead1.png
25+
0.75
26+
0.6
27+
Output.png
28+
OutputNearest.png
29+
)
30+
31+
if(ITK_WRAP_PYTHON)
32+
add_test(NAME ResampleSegmentedImageTestPython
33+
COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Code.py
34+
${CMAKE_CURRENT_BINARY_DIR}/2th_cthead1.png
35+
0.75
36+
0.6
37+
OutputPython.png
38+
OutputNearestPython.png
39+
)
40+
endif()

‎src/Core/ImageFunction/ResampleSegmentedImage/Code.cxx‎

Lines changed: 70 additions & 143 deletions
Original file line numberDiff line numberDiff line change
@@ -17,81 +17,41 @@
1717
*=========================================================================*/
1818
#include "itkImage.h"
1919
#include "itkImageFileReader.h"
20-
#include "itkIdentityTransform.h"
20+
#include "itkImageFileWriter.h"
2121
#include "itkLabelImageGaussianInterpolateImageFunction.h"
2222
#include "itkNearestNeighborInterpolateImageFunction.h"
2323
#include "itkResampleImageFilter.h"
2424

25-
#include "itkCustomColormapFunction.h"
26-
#include "itkScalarToRGBColormapImageFilter.h"
27-
#include "itkRGBPixel.h"
28-
#include "itkMersenneTwisterRandomVariateGenerator.h"
29-
30-
#ifdef ENABLE_QUICKVIEW
31-
# include "QuickView.h"
32-
#endif
33-
34-
#include <iostream>
35-
36-
using InputPixelType = unsigned short;
37-
using ImageType = itk::Image<InputPixelType, 2>;
38-
using RGBPixelType = itk::RGBPixel<unsigned char>;
39-
using RGBImageType = itk::Image<RGBPixelType, 2>;
40-
41-
using ColormapType = itk::Function::CustomColormapFunction<ImageType::PixelType, RGBImageType::PixelType>;
42-
43-
static void
44-
CreateRandomColormap(unsigned int size, ColormapType::Pointer colormap);
45-
4625
int
4726
main(int argc, char * argv[])
4827
{
49-
if (argc != 3)
28+
if (argc != 6)
5029
{
51-
std::cerr << "Usage: " << argv[0] << " inputImageFile pixelSize" << std::endl;
30+
std::cerr << "Usage: " << argv[0]
31+
<< " inputImageFile spacingFraction sigmaFraction outputImageFileLabelImageInterpolator "
32+
"outputImageFileNearestNeighborInterpolator"
33+
<< std::endl;
5234

5335
return EXIT_FAILURE;
5436
}
55-
double spacing = std::stod(argv[2]);
56-
57-
using ReaderType = itk::ImageFileReader<ImageType>;
58-
59-
// Identity transform.
60-
// We don't want any transform on our image except rescaling which is not
61-
// specified by a transform but by the input/output spacing as we will see
62-
// later.
63-
// So no transform will be specified.
64-
using TransformType = itk::IdentityTransform<double, 2>;
37+
const char * const inputImageFile = argv[1];
38+
const double spacingFraction = std::stod(argv[2]);
39+
const double sigmaFraction = std::stod(argv[3]);
40+
const char * const outputImageFileLabelImageInterpolator = argv[4];
41+
const char * const outputImageFileNearestNeighborInterpolator = argv[5];
6542

66-
using GaussianInterpolatorType = itk::LabelImageGaussianInterpolateImageFunction<ImageType, double>;
67-
using NearestNeighborInterpolatorType = itk::NearestNeighborInterpolateImageFunction<ImageType, double>;
43+
constexpr unsigned int Dimension = 2;
44+
using PixelType = unsigned char;
45+
using ImageType = itk::Image<PixelType, Dimension>;
6846

69-
using ResampleFilterType = itk::ResampleImageFilter<ImageType, ImageType>;
70-
71-
// Prepare the reader and update it right away to know the sizes beforehand.
47+
using ReaderType = itk::ImageFileReader<ImageType>;
7248
ReaderType::Pointer reader = ReaderType::New();
73-
reader->SetFileName(argv[1]);
49+
reader->SetFileName(inputImageFile);
7450
reader->Update();
7551

76-
// Instantiate the transform and specify it should be the identity transform.
77-
TransformType::Pointer transform = TransformType::New();
78-
transform->SetIdentity();
79-
80-
// Instantiate the interpolators
81-
GaussianInterpolatorType::Pointer gaussianInterpolator = GaussianInterpolatorType::New();
82-
gaussianInterpolator->SetSigma(1.0);
83-
gaussianInterpolator->SetAlpha(3.0);
84-
85-
NearestNeighborInterpolatorType::Pointer nearestNeighborInterpolator = NearestNeighborInterpolatorType::New();
86-
87-
// Instantiate the resamplers. Wire in the transforms and the interpolators.
88-
ResampleFilterType::Pointer resizeFilter1 = ResampleFilterType::New();
89-
resizeFilter1->SetTransform(transform);
90-
resizeFilter1->SetInterpolator(gaussianInterpolator);
91-
92-
ResampleFilterType::Pointer resizeFilter2 = ResampleFilterType::New();
93-
resizeFilter2->SetTransform(transform);
94-
resizeFilter2->SetInterpolator(nearestNeighborInterpolator);
52+
using ResampleFilterType = itk::ResampleImageFilter<ImageType, ImageType>;
53+
ResampleFilterType::Pointer resizeFilter = ResampleFilterType::New();
54+
resizeFilter->SetInput(reader->GetOutput());
9555

9656
// Compute and set the output size
9757
//
@@ -109,96 +69,63 @@ main(int argc, char * argv[])
10969
// So either we specify new height and width and compute new spacings
11070
// or we specify new spacing and compute new height and width
11171
// and computations that follows need to be modified a little (as it is
112-
// done at step 2 there:
113-
// http://itk.org/Wiki/ITK/Examples/DICOM/ResampleDICOM)
114-
//
115-
116-
// Fetch original image spacing.
117-
const ImageType::SpacingType & inputSpacing = reader->GetOutput()->GetSpacing();
118-
119-
double outputSpacing[2];
120-
outputSpacing[0] = spacing;
121-
outputSpacing[1] = spacing;
122-
123-
// Fetch original image size
124-
const ImageType::RegionType & inputRegion = reader->GetOutput()->GetLargestPossibleRegion();
125-
const ImageType::SizeType & inputSize = inputRegion.GetSize();
126-
unsigned int oldWidth = inputSize[0];
127-
unsigned int oldHeight = inputSize[1];
128-
129-
unsigned int newWidth = (double)oldWidth * inputSpacing[0] / spacing;
130-
unsigned int newHeight = (double)oldHeight * inputSpacing[1] / spacing;
131-
132-
// Set the output spacing as specified on the command line
133-
resizeFilter1->SetOutputSpacing(outputSpacing);
134-
resizeFilter2->SetOutputSpacing(outputSpacing);
135-
136-
// Set the computed size
137-
itk::Size<2> outputSize = { { newWidth, newHeight } };
138-
resizeFilter1->SetSize(outputSize);
139-
resizeFilter2->SetSize(outputSize);
14072

141-
// Specify the input for the resamplers
142-
resizeFilter1->SetInput(reader->GetOutput());
143-
resizeFilter2->SetInput(reader->GetOutput());
144-
145-
// Display the images
146-
using ColormapFilterType = itk::ScalarToRGBColormapImageFilter<ImageType, RGBImageType>;
147-
ColormapFilterType::Pointer colormapFilter1 = ColormapFilterType::New();
148-
149-
ColormapType::Pointer colormap = ColormapType::New();
150-
CreateRandomColormap(4096, colormap);
151-
152-
colormapFilter1->SetInput(reader->GetOutput());
153-
colormapFilter1->SetColormap(colormap);
154-
155-
ColormapFilterType::Pointer colormapFilter2 = ColormapFilterType::New();
156-
colormapFilter2->SetInput(resizeFilter1->GetOutput());
157-
colormapFilter2->SetColormap(colormap);
158-
159-
ColormapFilterType::Pointer colormapFilter3 = ColormapFilterType::New();
160-
colormapFilter3->SetInput(resizeFilter2->GetOutput());
161-
colormapFilter3->SetColormap(colormap);
73+
const ImageType::SpacingType inputSpacing{ reader->GetOutput()->GetSpacing() };
74+
ImageType::SpacingType outputSpacing;
75+
for (unsigned int dim = 0; dim < Dimension; ++dim)
76+
{
77+
outputSpacing[dim] = inputSpacing[dim] * spacingFraction;
78+
}
79+
resizeFilter->SetOutputSpacing(outputSpacing);
16280

163-
#ifdef ENABLE_QUICKVIEW
164-
std::stringstream desc;
165-
desc << itksys::SystemTools::GetFilenameName(argv[1]) << ": " << oldWidth << ", " << oldHeight;
166-
QuickView viewer;
167-
viewer.AddRGBImage(colormapFilter1->GetOutput(), true, desc.str());
81+
const ImageType::RegionType inputRegion{ reader->GetOutput()->GetLargestPossibleRegion() };
82+
const ImageType::SizeType inputSize{ inputRegion.GetSize() };
83+
ImageType::SizeType outputSize;
84+
for (unsigned int dim = 0; dim < Dimension; ++dim)
85+
{
86+
outputSize[dim] = inputSize[dim] * inputSpacing[dim] / spacingFraction;
87+
}
88+
resizeFilter->SetSize(outputSize);
16889

169-
std::stringstream desc2;
170-
desc2 << "Gaussian Interpolation: " << newWidth << ", " << newHeight;
171-
viewer.AddRGBImage(colormapFilter2->GetOutput(), true, desc2.str());
90+
using GaussianInterpolatorType = itk::LabelImageGaussianInterpolateImageFunction<ImageType, double>;
91+
GaussianInterpolatorType::Pointer gaussianInterpolator = GaussianInterpolatorType::New();
92+
GaussianInterpolatorType::ArrayType sigma;
93+
for (unsigned int dim = 0; dim < Dimension; ++dim)
94+
{
95+
sigma[dim] = outputSpacing[dim] * sigmaFraction;
96+
}
97+
gaussianInterpolator->SetSigma(sigma);
98+
gaussianInterpolator->SetAlpha(3.0);
99+
resizeFilter->SetInterpolator(gaussianInterpolator);
172100

173-
std::stringstream desc3;
174-
desc3 << "Nearest Neighbor Interpolation: " << newWidth << ", " << newHeight;
175-
viewer.AddRGBImage(colormapFilter3->GetOutput(), true, desc3.str());
176-
viewer.Visualize();
177-
#endif
178-
return EXIT_SUCCESS;
179-
}
101+
using WriterType = itk::ImageFileWriter<ImageType>;
102+
WriterType::Pointer writer = WriterType::New();
103+
writer->SetInput(resizeFilter->GetOutput());
104+
writer->SetFileName(outputImageFileLabelImageInterpolator);
105+
try
106+
{
107+
writer->Update();
108+
}
109+
catch (itk::ExceptionObject & error)
110+
{
111+
std::cerr << "Error: " << error << std::endl;
112+
return EXIT_FAILURE;
113+
}
180114

181-
void
182-
CreateRandomColormap(unsigned int size, ColormapType::Pointer colormap)
183-
{
184-
#define LOW .3
185-
ColormapType::ChannelType redChannel;
186-
ColormapType::ChannelType greenChannel;
187-
ColormapType::ChannelType blueChannel;
188-
itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer random =
189-
itk::Statistics::MersenneTwisterRandomVariateGenerator::New();
115+
using NearestNeighborInterpolatorType = itk::NearestNeighborInterpolateImageFunction<ImageType, double>;
116+
NearestNeighborInterpolatorType::Pointer nearestNeighborInterpolator = NearestNeighborInterpolatorType::New();
117+
resizeFilter->SetInterpolator(nearestNeighborInterpolator);
190118

191-
random->SetSeed(8775070);
192-
redChannel.push_back(LOW);
193-
greenChannel.push_back(LOW);
194-
blueChannel.push_back(LOW);
195-
for (unsigned int i = 1; i < size; ++i)
119+
writer->SetFileName(outputImageFileNearestNeighborInterpolator);
120+
try
196121
{
197-
redChannel.push_back(static_cast<ColormapType::RealType>(random->GetUniformVariate(LOW, 1.0)));
198-
greenChannel.push_back(static_cast<ColormapType::RealType>(random->GetUniformVariate(LOW, 1.0)));
199-
blueChannel.push_back(static_cast<ColormapType::RealType>(random->GetUniformVariate(LOW, 1.0)));
122+
writer->Update();
200123
}
201-
colormap->SetRedChannel(redChannel);
202-
colormap->SetGreenChannel(greenChannel);
203-
colormap->SetBlueChannel(blueChannel);
124+
catch (itk::ExceptionObject & error)
125+
{
126+
std::cerr << "Error: " << error << std::endl;
127+
return EXIT_FAILURE;
128+
}
129+
130+
return EXIT_SUCCESS;
204131
}

‎src/Core/ImageFunction/ResampleSegmentedImage/Code.py‎

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
#!/usr/bin/env python
2+
3+
# Copyright NumFOCUS
4+
#
5+
# Licensed under the Apache License, Version 2.0 (the "License");
6+
# you may not use this file except in compliance with the License.
7+
# You may obtain a copy of the License at
8+
#
9+
# http://www.apache.org/licenses/LICENSE-2.0.txt
10+
#
11+
# Unless required by applicable law or agreed to in writing, software
12+
# distributed under the License is distributed on an "AS IS" BASIS,
13+
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
# See the License for the specific language governing permissions and
15+
# limitations under the License.
16+
17+
import itk
18+
import sys
19+
20+
if len(sys.argv) != 6:
21+
print(f"Usage: {sys.argv[0]} input_image_file spacing_fraction sigma_fraction output_image_file_label_image_interpolator output_image_file_nearest_neighbor_interpolator")
22+
sys.exit(1)
23+
24+
input_image_file = sys.argv[1]
25+
spacing_fraction = float(sys.argv[2])
26+
sigma_fraction = float(sys.argv[3])
27+
output_image_file_label_image_interpolator = sys.argv[4]
28+
output_image_file_nearest_neighbor_interpolator = sys.argv[5]
29+
30+
input_image = itk.imread(input_image_file)
31+
32+
resize_filter = itk.ResampleImageFilter.New(input_image)
33+
34+
input_spacing = itk.spacing(input_image)
35+
output_spacing = [s * spacing_fraction for s in input_spacing]
36+
resize_filter.SetOutputSpacing(output_spacing)
37+
38+
input_size = itk.size(input_image)
39+
output_size = [int(s * input_spacing[dim] / spacing_fraction) for dim, s in
40+
enumerate(input_size)]
41+
resize_filter.SetSize(output_size)
42+
43+
gaussian_interpolator = itk.LabelImageGaussianInterpolateImageFunction.New(input_image)
44+
sigma = [s * sigma_fraction for s in output_spacing]
45+
gaussian_interpolator.SetSigma(sigma)
46+
gaussian_interpolator.SetAlpha(3.0)
47+
resize_filter.SetInterpolator(gaussian_interpolator)
48+
49+
itk.imwrite(resize_filter, output_image_file_label_image_interpolator)
50+
51+
nearest_neighbor_interpolator = itk.NearestNeighborInterpolateImageFunction.New(input_image)
52+
resize_filter.SetInterpolator(nearest_neighbor_interpolator)
53+
54+
itk.imwrite(resize_filter, output_image_file_nearest_neighbor_interpolator)

0 commit comments

Comments
 (0)