From 0ca3a0f02f145bcef32706c1e77ed8e2b79d125d Mon Sep 17 00:00:00 2001 From: YuanhaoGong Date: Sat, 24 Oct 2015 14:56:03 +0800 Subject: [PATCH] Add ITK version --- DM.h | 15 ++++ ITK/CMakeLists.txt | 8 +++ ITK/itkCurvatureFilter.h | 150 +++++++++++++++++++++++++++++++++++++++ ITK/main.cxx | 99 ++++++++++++++++++++++++++ ReadMe.txt | 2 + main.cpp | 19 ++++- main_MultiScale.cpp | 15 ++++ 7 files changed, 306 insertions(+), 2 deletions(-) create mode 100644 ITK/CMakeLists.txt create mode 100644 ITK/itkCurvatureFilter.h create mode 100644 ITK/main.cxx diff --git a/DM.h b/DM.h index fd0c472..8c98668 100644 --- a/DM.h +++ b/DM.h @@ -1,3 +1,18 @@ +/*========================================================================= + * + * Curvature Filter + * + ************************************************************************** + + @phdthesis{gong:phd, + title={Spectrally regularized surfaces}, + author={Gong, Yuanhao}, + year={2015}, + school={ETH Zurich, Nr. 22616}, + note={http://dx.doi.org/10.3929/ethz-a-010438292}} + + *=========================================================================*/ + //Dual Mesh sctructure class DM { diff --git a/ITK/CMakeLists.txt b/ITK/CMakeLists.txt new file mode 100644 index 0000000..6ea6465 --- /dev/null +++ b/ITK/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 2.8) +set(CMAKE_CONFIGURATION_TYPES Release) +project( CF ) +find_package( ITK REQUIRED ) +include( ${ITK_USE_FILE} ) +add_executable( CF main.cxx ) +target_link_libraries( CF ${ITK_LIBRARIES} ) +set(CMAKE_BUILD_TYPE Release) diff --git a/ITK/itkCurvatureFilter.h b/ITK/itkCurvatureFilter.h new file mode 100644 index 0000000..ca0202c --- /dev/null +++ b/ITK/itkCurvatureFilter.h @@ -0,0 +1,150 @@ +/*========================================================================= + * + * Curvature Filter in ITK, please cite following work in you paper! + * + ************************************************************************** + + @phdthesis{gong:phd, + title={Spectrally regularized surfaces}, + author={Gong, Yuanhao}, + year={2015}, + school={ETH Zurich, Nr. 22616}, + note={http://dx.doi.org/10.3929/ethz-a-010438292}} + + *=========================================================================*/ +//estimate the d_m by different curvature filters +float GC(float & cur, float & pre, float & rigUp, float & right, float & rigDn, float & next, float & lefDn, float & left, float & lefUp); +float MC(float & cur, float & pre, float & rigUp, float & right, float & rigDn, float & next, float & lefDn, float & left, float & lefUp); +float TV(float & cur, float & pre, float & rigUp, float & right, float & rigDn, float & next, float & lefDn, float & left, float & lefUp); +float (*projection)(float & cur, float & pre, float & rigUp, float & right, float & rigDn, float & next, float & lefDn, float & left, float & lefUp); +void CurvatureFilter(itk::Image< float, 2 >::Pointer image, int FilterType, int IterationNumber) +{ + + if (FilterType == 0) {projection = TV; std::cout<<"TV filter: with "<::SizeType size = image->GetLargestPossibleRegion().GetSize(); + itk::Image< float, 2 >::IndexType index_cur, index_prev, index_rigUp, index_right, index_rigDn, index_next, index_leftDn, index_left, index_leftUp; + + clock_t Tstart, Tend; + Tstart = clock(); + + for(int it = 0; it < IterationNumber; ++it) + { + //domain decomposition, four sets + for (int i = 1; i < size[0]; ++i, ++i) + for (int j = 1; j < size[1]; ++j, ++j) + { + index_leftUp[0] = j-1, index_leftUp[1] = i-1, index_prev[0] = j, index_prev[1] = i-1, index_rigUp[0]=j+1, index_rigUp[1] = i-1; + index_left[0]=j-1, index_left[1]=i, index_cur[0] = j, index_cur[1] = i, index_right[0]=j+1, index_right[1]=i; + index_leftDn[0]=j-1, index_leftDn[1]=i+1, index_next[0]=j, index_next[1] = i+1, index_rigDn[0]=j+1, index_rigDn[1] = i+1; + + cur = image->GetPixel(index_cur); + d = (*projection)(cur, image->GetPixel(index_prev), image->GetPixel(index_rigUp), image->GetPixel(index_right), image->GetPixel(index_rigDn), image->GetPixel(index_next), image->GetPixel(index_leftDn), image->GetPixel(index_left), image->GetPixel(index_leftUp)); + image->SetPixel(index_cur, cur + d); + } + for (int i = 2; i < size[0]; ++i, ++i) + for (int j = 2; j < size[1]; ++j, ++j) + { + index_leftUp[0] = j-1, index_leftUp[1] = i-1, index_prev[0] = j, index_prev[1] = i-1, index_rigUp[0]=j+1, index_rigUp[1] = i-1; + index_left[0]=j-1, index_left[1]=i, index_cur[0] = j, index_cur[1] = i, index_right[0]=j+1, index_right[1]=i; + index_leftDn[0]=j-1, index_leftDn[1]=i+1, index_next[0]=j, index_next[1] = i+1, index_rigDn[0]=j+1, index_rigDn[1] = i+1; + + cur = image->GetPixel(index_cur); + d = (*projection)(cur, image->GetPixel(index_prev), image->GetPixel(index_rigUp), image->GetPixel(index_right), image->GetPixel(index_rigDn), image->GetPixel(index_next), image->GetPixel(index_leftDn), image->GetPixel(index_left), image->GetPixel(index_leftUp)); + image->SetPixel(index_cur, cur + d); + } + for (int i = 1; i < size[0]; ++i, ++i) + for (int j = 2; j < size[1]; ++j, ++j) + { + index_leftUp[0] = j-1, index_leftUp[1] = i-1, index_prev[0] = j, index_prev[1] = i-1, index_rigUp[0]=j+1, index_rigUp[1] = i-1; + index_left[0]=j-1, index_left[1]=i, index_cur[0] = j, index_cur[1] = i, index_right[0]=j+1, index_right[1]=i; + index_leftDn[0]=j-1, index_leftDn[1]=i+1, index_next[0]=j, index_next[1] = i+1, index_rigDn[0]=j+1, index_rigDn[1] = i+1; + + cur = image->GetPixel(index_cur); + d = (*projection)(cur, image->GetPixel(index_prev), image->GetPixel(index_rigUp), image->GetPixel(index_right), image->GetPixel(index_rigDn), image->GetPixel(index_next), image->GetPixel(index_leftDn), image->GetPixel(index_left), image->GetPixel(index_leftUp)); + image->SetPixel(index_cur, cur + d); + } + for (int i = 2; i < size[0]; ++i, ++i) + for (int j = 1; j < size[1]; ++j, ++j) + { + index_leftUp[0] = j-1, index_leftUp[1] = i-1, index_prev[0] = j, index_prev[1] = i-1, index_rigUp[0]=j+1, index_rigUp[1] = i-1; + index_left[0]=j-1, index_left[1]=i, index_cur[0] = j, index_cur[1] = i, index_right[0]=j+1, index_right[1]=i; + index_leftDn[0]=j-1, index_leftDn[1]=i+1, index_next[0]=j, index_next[1] = i+1, index_rigDn[0]=j+1, index_rigDn[1] = i+1; + + cur = image->GetPixel(index_cur); + d = (*projection)(cur, image->GetPixel(index_prev), image->GetPixel(index_rigUp), image->GetPixel(index_right), image->GetPixel(index_rigDn), image->GetPixel(index_next), image->GetPixel(index_leftDn), image->GetPixel(index_left), image->GetPixel(index_leftUp)); + image->SetPixel(index_cur, cur + d); + } + } + + Tend = clock() - Tstart; + double time = double(Tend)/(CLOCKS_PER_SEC/1000.0); + std::cout<<" running time is "< +#include "itkCurvatureFilter.h" + +int FType = 2; + +int main( int argc, char * argv[] ) +{ + if( argc < 5 ) + { + std::cerr << "Usage: " << std::endl; + std::cerr << argv[0] << " inputImageFile outputImageFile filterType numberOfIterations" << std::endl; + return EXIT_FAILURE; + } + + const unsigned int numberOfIterations = atoi( argv[4] ); + const char * filterType = argv[3]; + + if (*filterType == 't') {FType = 0; } + if (*filterType == 'm') {FType = 1; } + if (*filterType == 'g') {FType = 2; } + + // read image + typedef float PixelType; + typedef itk::Image< PixelType, 2 > ImageType; + typedef itk::ImageFileReader< ImageType > ReaderType; + + typedef itk::ConstNeighborhoodIterator< ImageType > NeighborhoodIteratorType; + typedef itk::ImageRegionIterator< ImageType> IteratorType; + + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName( argv[1] ); + try + { + reader->Update(); + } + catch ( itk::ExceptionObject &err) + { + std::cout << "ExceptionObject caught !" << std::endl; + std::cout << err << std::endl; + return -1; + } + ImageType::Pointer image = reader->GetOutput(); + + /*********** curvature filter ******************/ + CurvatureFilter(image, FType, numberOfIterations); + /*********** curvature filter ******************/ + + // save the result + typedef unsigned char WritePixelType; + typedef itk::Image< WritePixelType, 2 > WriteImageType; + typedef itk::ImageFileWriter< WriteImageType > WriterType; + + typedef itk::RescaleIntensityImageFilter< + ImageType, WriteImageType > RescaleFilterType; + + RescaleFilterType::Pointer rescaler = RescaleFilterType::New(); + + rescaler->SetOutputMinimum( 0 ); + rescaler->SetOutputMaximum( 255 ); + rescaler->SetInput(image); + + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName( argv[2] ); + writer->SetInput(rescaler->GetOutput()); + try + { + writer->Update(); + } + catch ( itk::ExceptionObject &err) + { + std::cout << "ExceptionObject caught !" << std::endl; + std::cout << err << std::endl; + return -1; + } + + return EXIT_SUCCESS; +} + diff --git a/ReadMe.txt b/ReadMe.txt index 88320d5..7a84b22 100644 --- a/ReadMe.txt +++ b/ReadMe.txt @@ -3,12 +3,14 @@ This code was developed by Yuanhao Gong during his PhD at MOSAIC Group. Please cite Yuanhao's PhD thesis if you use this code in your work. Thank you! ============================================================= + @phdthesis{gong:phd, title={Spectrally regularized surfaces}, author={Gong, Yuanhao}, year={2015}, school={ETH Zurich, Nr. 22616}, note={http://dx.doi.org/10.3929/ethz-a-010438292}} + ============================================================= FAQ: diff --git a/main.cpp b/main.cpp index 8a35bc5..45220fb 100644 --- a/main.cpp +++ b/main.cpp @@ -1,3 +1,18 @@ +/*========================================================================= + * + * Curvature Filter + * + ************************************************************************** + + @phdthesis{gong:phd, + title={Spectrally regularized surfaces}, + author={Gong, Yuanhao}, + year={2015}, + school={ETH Zurich, Nr. 22616}, + note={http://dx.doi.org/10.3929/ethz-a-010438292}} + + *=========================================================================*/ + #include #include #include @@ -71,8 +86,8 @@ int main(int argc, char** argv) if (argc==6) { - //filter solver for the variational models - DualMesh.read(argv[1]); + //filter solver for the variational models + DualMesh.read(argv[1]); DualMesh.Solver(Type, mytime, ItNum, lambda, DataFitOrder); cout<<"runtime is "< #include #include