Add ITK version

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
 {

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)

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 "<<IterationNumber<<" Iterations";}
+ if (FilterType == 1) {projection = MC; std::cout<<"MC filter: with "<<IterationNumber<<" Iterations";}
+ if (FilterType == 2) {projection = GC; std::cout<<"GC filter: with "<<IterationNumber<<" Iterations";}
+ float d, cur; 
+ //image size and neighbor index
+ itk::Image< float, 2 >::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 "<<time<<" milliseconds"<<std::endl;
+}
+
+// *************************** functions ********************************************
+//estimate d_m (minimal projection signed distance)
+float GC(float & cur, float & prev, float & rigUp, float & right, float & rigDn, float & down, float & lefDn, float & left, float& lefUp)
+{
+ float d_m = (prev + down)/2 - cur;
+ float tmp = (left, right)/2 - cur;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+ tmp = (lefUp, rigDn)/2 - cur;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+ tmp = (lefDn, rigUp)/2 - cur;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+
+ //hybrid
+ tmp = (left + lefUp + prev)/3 - cur;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+ tmp = (prev + rigUp + right)/3 - cur;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+ tmp = (right + rigDn + down)/3 - cur;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp; 
+ tmp = (down + lefDn + left)/3 - cur;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+
+ return d_m;
+}
+float MC(float & cur, float & prev, float & rigUp, float & right, float & rigDn, float & down, float & lefDn, float & left, float& lefUp)
+{
+ float TM = 8*cur;
+
+ float d_m = (prev + down)*2.5f - rigUp - rigDn + 5*right - TM;
+ float tmp = (prev + down)*2.5f - lefUp - lefDn + 5*left - TM;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+ tmp = (left + right)*2.5f - lefUp - rigUp + 5*prev - TM;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+ tmp = (left + right)*2.5f - lefDn - rigDn + 5*down - TM;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+
+ return d_m/8;
+}
+float TV(float & cur, float & prev, float & rigUp, float & right, float & rigDn, float & down, float & lefDn, float & left, float& lefUp)
+{
+ float TM = 5*cur;
+
+ float d_m = prev + down + lefUp + left + lefDn - TM;
+ float tmp = prev + down + rigUp + rigDn + right - TM;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+ tmp = left + lefUp + prev + rigUp + right - TM;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+ tmp = left + lefDn + down + rigDn + right - TM;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+
+ //diag
+ tmp = lefUp + prev + rigUp + left + lefDn - TM;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+ tmp = lefUp + prev + rigUp + right + rigDn - TM;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+ tmp = lefDn + down + rigDn + left + lefUp - TM;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+ tmp = lefDn + down + rigDn + right + rigUp - TM;
+ if (fabs(tmp) < fabs(d_m)) d_m = tmp;
+
+ return d_m/5;
+}
+

ITK/main.cxx

@@ -0,0 +1,99 @@
+/*=========================================================================
+ *
+ * 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}}
+
+ *=========================================================================*/
+
+#include "itkImageFileReader.h"
+#include "itkImageFileWriter.h"
+#include "itkRescaleIntensityImageFilter.h"
+
+#include "itkConstNeighborhoodIterator.h"
+#include "itkImageRegionIterator.h"
+#include <time.h>
+#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;
+}
+

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:

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 <opencv2/core/core.hpp>
 #include <opencv2/highgui/highgui.hpp>
 #include <opencv2/imgproc/imgproc.hpp>
@@ -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 "<<mytime<<" milliseconds."<<endl;
  DualMesh.write("CF_Solver.png");

main_MultiScale.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 <opencv2/core/core.hpp>
 #include <opencv2/highgui/highgui.hpp>
 #include <opencv2/imgproc/imgproc.hpp>