1、代码
1 #include "itkImage.h"
2 #include "itkImageFileReader.h"
3 #include "itkSignedMaurerDistanceMapImageFilter.h"
4
5 #include "itksys/SystemTools.hxx"
6 #include <sstream>
7 //#define ENABLE_QUICKVIEW
8 #ifdef ENABLE_QUICKVIEW
9 #include "QuickView.h"
10 #endif
11
12 using UnsignedCharImageType = itk::Image<unsigned char, 2>;
13 using FloatImageType = itk::Image<float, 2>;
14
15 static void
16 CreateImage(UnsignedCharImageType::Pointer image);
17
18 int
19 main(int argc, char * argv[])
20 {
21 UnsignedCharImageType::Pointer image = UnsignedCharImageType::New();
22 if (argc < 2)
23 {
24 CreateImage(image);
25 }
26 else
27 {
28 using ReaderType = itk::ImageFileReader<UnsignedCharImageType>;
29 ReaderType::Pointer reader = ReaderType::New();
30 reader->SetFileName(argv[1]);
31 reader->Update();
32 image = reader->GetOutput();
33 }
34
35 using SignedMaurerDistanceMapImageFilterType =
36 itk::SignedMaurerDistanceMapImageFilter<UnsignedCharImageType, FloatImageType>;
37 SignedMaurerDistanceMapImageFilterType::Pointer distanceMapImageFilter =
38 SignedMaurerDistanceMapImageFilterType::New();
39 distanceMapImageFilter->SetInput(image);
40 distanceMapImageFilter->Update();
41 auto outputImage = distanceMapImageFilter->GetOutput();
42
43 for (unsigned int i = 40; i < 50; ++i)
44 {
45 for (int j = 40; j < 50; j++)
46 {
47 itk::Index<2> pixel{i,j};
48 //pixel.Fill(i);
49 std::cout << outputImage->GetPixel(pixel) << "\t";
50 }
51 std::cout << std::endl;
52 }
53
54 #ifdef ENABLE_QUICKVIEW
55 QuickView viewer;
56 viewer.AddImage(
57 image.GetPointer(), true, argc > 1 ? itksys::SystemTools::GetFilenameName(argv[1]) : "Generated image");
58
59 std::stringstream desc;
60 desc << "Signed Maurer Distance";
61 viewer.AddImage(distanceMapImageFilter->GetOutput(), true, desc.str());
62
63 viewer.Visualize();
64 #endif
65
66 return EXIT_SUCCESS;
67 }
68
69
70 void
71 CreateImage(UnsignedCharImageType::Pointer image)
72 {
73 // Create an image
74 itk::Index<2> start;
75 start.Fill(0);
76
77 itk::Size<2> size;
78 size.Fill(100);
79
80 itk::ImageRegion<2> region(start, size);
81 image->SetRegions(region);
82 image->Allocate();
83 image->FillBuffer(0);
84
85 // Create a line of white pixels
86 for (unsigned int i = 40; i < 60; ++i)
87 {
88 itk::Index<2> pixel;
89 pixel.Fill(i);
90 image->SetPixel(pixel, 255);
91 }
92 }
2、CMakeLists.txt
关于 CMakeLists.txt 可以修改只保留itk,下面是因为我测试过itk和vtk一起使用才写成下面的样子的。
1 cmake_minimum_required(VERSION 2.8 FATAL_ERROR)
2 project(testfilter)
3
4 set(VTK_DIR D:/ProgramFiles/vtk-9.0/lib/cmake/vtk-9.0)
5 find_package(VTK COMPONENTS
6 vtkCommonCore
7 vtkCommonDataModel
8 vtkFiltersSources
9 vtkIOImage
10 vtkImagingStencil
11 vtkFiltersCore
12 vtkRenderingCore
13 vtkFiltersHybrid
14 vtkInteractionStyle
15 vtkRenderingFreeType
16 vtkRenderingOpenGL2
17 vtkRenderingVolumeOpenGL2
18 vtkImagingColor
19 QUIET
20 )
21
22 if (NOT VTK_FOUND)
23 message("Skipping testfilter: ${VTK_NOT_FOUND_MESSAGE}")
24 return()
25 endif()
26 message (STATUS "VTK_VERSION: ${VTK_VERSION}")
27
28 if (VTK_VERSION VERSION_LESS "8.90.0")
29 # old system
30 include(${VTK_USE_FILE})
31 endif()
32
33 set(ITK_DIR D:/ProgramFiles/ITK-5.2/lib/cmake/ITK-5.2)
34
35 FIND_PACKAGE(ITK REQUIRED)
36 INCLUDE(${ITK_USE_FILE})
37
38 file(GLOB_RECURSE mains ${CMAKE_CURRENT_SOURCE_DIR} *.cpp)
39 message(STATUS ${mains})
40 foreach(mainfile IN LISTS mains)
41 # Get file name without directory
42 get_filename_component(mainname ${mainfile} NAME_WE)
43 add_executable(${mainname} ${mainfile})
44 target_link_libraries (${mainname} ${VTK_LIBRARIES} ${ITK_LIBRARIES})
45 # vtk_module_autoinit is needed
46 if (NOT VTK_VERSION VERSION_LESS "8.90.0")
47 # vtk_module_autoinit is needed
48 vtk_module_autoinit(
49 TARGETS ${mainname}
50 MODULES ${VTK_LIBRARIES}
51 )
52 endif()
53
54 endforeach()
3、测试输出
下面的输出是我们想要的结果
1 -0 1 1.41421 2.23607 2.82843 3.60555 4.24264 5 5.65685 6.40312
2 1 -0 1 1.41421 2.23607 2.82843 3.60555 4.24264 5 5.65685
3 1.41421 1 -0 1 1.41421 2.23607 2.82843 3.60555 4.24264 5
4 2.23607 1.41421 1 -0 1 1.41421 2.23607 2.82843 3.60555 4.24264
5 2.82843 2.23607 1.41421 1 -0 1 1.41421 2.23607 2.82843 3.60555
6 3.60555 2.82843 2.23607 1.41421 1 -0 1 1.41421 2.23607 2.82843
7 4.24264 3.60555 2.82843 2.23607 1.41421 1 -0 1 1.41421 2.23607
8 5 4.24264 3.60555 2.82843 2.23607 1.41421 1 -0 1 1.41421
9 5.65685 5 4.24264 3.60555 2.82843 2.23607 1.41421 1 -0 1
10 6.40312 5.65685 5 4.24264 3.60555 2.82843 2.23607 1.41421 1 -0
该结果与scipy中ndimage中的distance_transform_edt函数作用相同
ndimage.distance_transform_edt(seed == 0)