changeset 25:6c6003bcad16

Cleanup, work on VTK code
author Jordi Gutiérrez Hermoso <jordigh@gmail.com>
date Mon, 01 Feb 2010 17:42:16 -0600
parents ad9e3d28ce9b
children 093a9e7e26df
files Jamroot Makefile linalg.cpp main.cpp vtk/CMakeLists.txt vtk/terrain_test.cpp
diffstat 6 files changed, 238 insertions(+), 164 deletions(-) [+]
line wrap: on
line diff
deleted file mode 100644
--- a/Jamroot
+++ /dev/null
@@ -1,5 +0,0 @@
-lib gsl : : <name>gsl ;
-lib gslcblas : : <name>gslcblas ;
-exe rbf : bvp.cpp ddm.cpp diff_op.cpp error.cpp func.cpp interpolator.cpp linalg.cpp main.cpp rbf.cpp utils.cpp gsl gslcblas 
-	; 
-
--- a/Makefile
+++ b/Makefile
@@ -45,4 +45,7 @@
 
 cleandata:
 	rm -f data/*.matrix data/*.map data/*.vector \
-	 	results/*.matrix results/*.map results/*.png results/*.vector
\ No newline at end of file
+	 	results/*.matrix results/*.map results/*.png results/*.vector
+# DO NOT DELETE
+
+
--- a/linalg.cpp
+++ b/linalg.cpp
@@ -62,21 +62,21 @@
       gsl_matrix_memcpy(A,M.A);
 
       if(LUfactored)
-	delete LUptr;
+        delete LUptr;
 
       LUfactored = M.LUfactored;
       if(LUfactored){
-	LUptr = new LUmatrix(*M.LUptr);
+        LUptr = new LUmatrix(*M.LUptr);
       }
       
       if(SVDfactored)
-	gsl_vector_free(SVD);
+        gsl_vector_free(SVD);
 
       SVDfactored = M.SVDfactored;
       if(SVDfactored){
-	condition_number = M.condition_number;
-	SVD = gsl_vector_alloc(M.SVD->size);
-	gsl_vector_memcpy(SVD,M.SVD);
+        condition_number = M.condition_number;
+        SVD = gsl_vector_alloc(M.SVD->size);
+        gsl_vector_memcpy(SVD,M.SVD);
       }       
     }
     return *this;
@@ -231,7 +231,7 @@
     if(!LUfactored){
       LUptr = new LUmatrix(A);
       gsl_linalg_LU_decomp(LUptr->matrix_ptr(), LUptr->perm_ptr(), 
-			   LUptr->sgn_ptr());
+                           LUptr->sgn_ptr());
       LUfactored = true;
     }
     return LUptr;    
@@ -261,7 +261,7 @@
     matrix Z(A->size2, A->size1);
     for(size_t i = 1; i <= A->size1; i++)
       for(size_t j = 1; j <= A->size2; j++)
-	Z(j,i) = gsl_matrix_get(A,i-1,j-1);
+        Z(j,i) = gsl_matrix_get(A,i-1,j-1);
     return Z;
   }
   
@@ -323,7 +323,7 @@
 
     SVDfactor();
     condition_number =  gsl_vector_get(SVD,0) 
-                        /gsl_vector_get(SVD,std::max(A->size1, A->size2)-1);
+      /gsl_vector_get(SVD,std::max(A->size1, A->size2)-1);
     return condition_number;
   }
 
@@ -609,15 +609,15 @@
 
   bool vector::operator<(const vector& w) const{
     if(x -> size < w.x -> size) //Smaller vectors go first in this order.
-       return true;
+      return true;
 
     for(size_t i = 0; i < x -> size; i++){
       double L = gsl_vector_get(x,i);
       double R = gsl_vector_get(w.x,i);
       if(L < R )
-	return true;
+        return true;
       if(L > R )
-	return false;
+        return false;
     }
    
     return false; //Then vectors are equal.
@@ -661,7 +661,7 @@
     x -> size = (a.end() - a.begin())/a.stride()+1;
     x -> data = gsl_matrix_ptr(A, a.begin(), j-1);
     x -> stride = a.stride()*(A -> tda); //Note that a longer step is
-					 //necessary here.
+    //necessary here.
     x -> block = 0;
     x -> owner = 0;
   }
@@ -733,23 +733,23 @@
     if(b < a){
       badArgument exc;
       exc.reason = 
-	"Invalid arguments to slice::set. "
-	"Right endpoint must be greater than left endpoint.";
+        "Invalid arguments to slice::set. "
+        "Right endpoint must be greater than left endpoint.";
       exc.line = __LINE__; exc.file = __FILE__;
       throw exc;
     }
     else if(k==0){
       badArgument exc;
       exc.reason = 
-	"Invalid arguments to slice::set. "
-	"Cannot take zero stride.";
+        "Invalid arguments to slice::set. "
+        "Cannot take zero stride.";
       exc.line = __LINE__; exc.file = __FILE__;
       throw exc;
     }
     else if( k > b-a){
       badArgument exc;
       exc.reason = "Invalid arguments to slice::set. " 
-	"Step size cannot be greater than range.";
+        "Step size cannot be greater than range.";
       exc.line = __LINE__; exc.file = __FILE__;
       throw exc;
     }
@@ -772,71 +772,71 @@
     for(size_t i = 1; i <= v.size(); i++){
       os << " " << std::setw(int(v.precision()+6) ) << v(i) ;
     }
-     return os;
-   }
+    return os;
+  }
 
-   vector operator>>(std::istream& is, vector& v){
-     using namespace std;
-     string s;
-     list<double> data;
-     bool colvector = true;
-     bool shouldbedone = false;
-     while(getline(is, s)){
-       s = utils::trim(s);
-       if(s[0] == '#' or s.size() == 0) //Blank line or comment character
-	 continue;
+  vector operator>>(std::istream& is, vector& v){
+    using namespace std;
+    string s;
+    list<double> data;
+    bool colvector = true;
+    bool shouldbedone = false;
+    while(getline(is, s)){
+      s = utils::trim(s);
+      if(s[0] == '#' or s.size() == 0) //Blank line or comment character
+        continue;
 
-       stringstream ss;
-       ss << s;
+      stringstream ss;
+      ss << s;
 
-       double x;
-       size_t i = 0;
-       while(ss >> x){
-	 if( (i > 0 and colvector == false) or (shouldbedone == true)){
-	   badArgument exc;
-	   exc.reason = "Cannot read vector: bad format in input";
-	   exc.file = __FILE__;
-	   exc.line = __LINE__;
-	   throw exc;
-	 }
-	 data.push_back(x);
-	 i++;
-       }
-       if(i > 1){
-	 colvector = false; //So it must be a row vector instead
-	 shouldbedone = true;
-       }
-     }
+      double x;
+      size_t i = 0;
+      while(ss >> x){
+        if( (i > 0 and colvector == false) or (shouldbedone == true)){
+          badArgument exc;
+          exc.reason = "Cannot read vector: bad format in input";
+          exc.file = __FILE__;
+          exc.line = __LINE__;
+          throw exc;
+        }
+        data.push_back(x);
+        i++;
+      }
+      if(i > 1){
+        colvector = false; //So it must be a row vector instead
+        shouldbedone = true;
+      }
+    }
 
-     if(data.size() == 0){
-       endOfFile exc;
-       exc.reason = "Cannot read empty vector from input.";
-       exc.file = __FILE__;
-       exc.line = __LINE__;
-       throw exc;
-     }
+    if(data.size() == 0){
+      endOfFile exc;
+      exc.reason = "Cannot read empty vector from input.";
+      exc.file = __FILE__;
+      exc.line = __LINE__;
+      throw exc;
+    }
 
-     vector w(data.size());
-     typedef list<double>::iterator LI;
-     size_t k = 1;
-     for(LI i = data.begin(); i != data.end(); i++){
-       w(k) = *i;
-       k++;
-     }
-     v = w;
-     return v;
-   }
+    vector w(data.size());
+    typedef list<double>::iterator LI;
+    size_t k = 1;
+    for(LI i = data.begin(); i != data.end(); i++){
+      w(k) = *i;
+      k++;
+    }
+    v = w;
+    return v;
+  }
 
-   std::ostream& operator<<(std::ostream& os, const matrix& A){
-     os.setf(std::ios::scientific);
-     os << std::setprecision(int(A.precision()) );
-     for(size_t i = 1; i <= A.rows(); i++){
-       for(size_t j = 1; j <= A.cols(); j++)
-	 os << " " << std::setw(int(A.precision()+6) ) << A(i,j) << " ";
-       os << std::endl;
-     }
-     return os;
-   }
+  std::ostream& operator<<(std::ostream& os, const matrix& A){
+    os.setf(std::ios::scientific);
+    os << std::setprecision(int(A.precision()) );
+    for(size_t i = 1; i <= A.rows(); i++){
+      for(size_t j = 1; j <= A.cols(); j++)
+        os << " " << std::setw(int(A.precision()+6) ) << A(i,j) << " ";
+      os << std::endl;
+    }
+    return os;
+  }
 
   matrix operator>>(std::istream& is, matrix& A){
     using namespace std;
@@ -851,66 +851,66 @@
       line = utils::trim(line);
       //Blank row or comment character.
       if(line[0] == '#' or line.length() == 0) 
-	continue;
+        continue;
          
       stringstream ss_line;      
       cols = 0;
       ss_line << line;
       while(ss_line >> token){
-	if(token[0] == '#'){
-	  break; //Rest of line is comment.
-	}
+        if(token[0] == '#'){
+          break; //Rest of line is comment.
+        }
 
-	//The following may fail on a C++ implementation that doesn't
-	//obey IEEE arithmetic (IEC 559). We could check for those,
-	//but do we really want to compile Octave on C++
-	//implementations that don't follow IEEE arithmetic?
-	else if(token == "NaN"){ 
-	  double x = std::numeric_limits<double>::quiet_NaN();
-	  data.push_back(x);
-	  cols++;
-	}
-	else if(token == "Inf"){
-	  double x = std::numeric_limits<double>::infinity();
-	  data.push_back(x);
-	  cols++;
-	}
-	else if(token == "-Inf"){
-	  double x = -std::numeric_limits<double>::infinity();
-	  data.push_back(x);
-	  cols++;
-	}
-	else if(token == ","){
-	  ss_line >> token;
-	}
+        //The following may fail on a C++ implementation that doesn't
+        //obey IEEE arithmetic (IEC 559). We could check for those,
+        //but do we really want to compile Octave on C++
+        //implementations that don't follow IEEE arithmetic?
+        else if(token == "NaN"){ 
+          double x = std::numeric_limits<double>::quiet_NaN();
+          data.push_back(x);
+          cols++;
+        }
+        else if(token == "Inf"){
+          double x = std::numeric_limits<double>::infinity();
+          data.push_back(x);
+          cols++;
+        }
+        else if(token == "-Inf"){
+          double x = -std::numeric_limits<double>::infinity();
+          data.push_back(x);
+          cols++;
+        }
+        else if(token == ","){
+          ss_line >> token;
+        }
 
-	//This also ignores commas and any other token. I think. If
-	//there's garbage in the token, I have to see what happens
-	//here. Do we also need to check for garbage?
-	else{ 
-	  double x;
-	  stringstream ss_token;
-	  ss_token << token;
-	  ss_token >> x;
-	  data.push_back(x);
-	  cols++;
-	}
+        //This also ignores commas and any other token. I think. If
+        //there's garbage in the token, I have to see what happens
+        //here. Do we also need to check for garbage?
+        else{ 
+          double x;
+          stringstream ss_token;
+          ss_token << token;
+          ss_token >> x;
+          data.push_back(x);
+          cols++;
+        }
 
       }
        
       //First row gives the number of columns, and all successive rows
       //must have the same number of elements. 
       if(!rowset){
-	rowset = true;
-	rowsize = cols;
+        rowset = true;
+        rowsize = cols;
       }
        
       if (cols != rowsize){
-	badArgument exc;
-	exc.reason = "Cannot read matrix: bad format in input";
-	exc.file = __FILE__;
-	exc.line = __LINE__;
-	throw exc;
+        badArgument exc;
+        exc.reason = "Cannot read matrix: bad format in input";
+        exc.file = __FILE__;
+        exc.line = __LINE__;
+        throw exc;
       }
       rows++;
     }
@@ -928,8 +928,8 @@
     LI k = data.begin();
     for(size_t i = 1; i <= rows; i++){
       for(size_t j = 1; j <= cols; j++){
-	M(i,j) = *k;
-	k++;
+        M(i,j) = *k;
+        k++;
       }
     }
 
--- a/main.cpp
+++ b/main.cpp
@@ -47,7 +47,8 @@
 {  
   gsl_set_error_handler(&error_handling::errorHandler);
 
-  try{
+  try
+  {
     using namespace std;
     using boost::shared_ptr;
    
@@ -83,7 +84,8 @@
     }
 
   }
-  catch(error exc){
+  catch(error exc)
+  {
     utils::show_exception(exc);
     return 1;
   }
new file mode 100644
--- /dev/null
+++ b/vtk/CMakeLists.txt
@@ -0,0 +1,16 @@
+project (vtk_test)
+
+#set(CMAKE_VERBOSE_MAKEFILE TRUE)
+set(CMAKE_CXX_COMPILER "g++-4.4")
+
+cmake_minimum_required(VERSION 2.6)
+
+include (${CMAKE_ROOT}/Modules/FindVTK.cmake)
+
+if (USE_VTK_FILE)
+ include(${USE_VTK_FILE})
+endif (USE_VTK_FILE)
+
+add_executable(terrain_test terrain_test.cpp)
+target_link_libraries(terrain_test vtkRendering vtkHybrid)
+set_target_properties(terrain_test PROPERTIES COMPILE_FLAGS "-std=c++0x")
--- a/vtk/terrain_test.cpp
+++ b/vtk/terrain_test.cpp
@@ -19,19 +19,47 @@
 
 #include <vtkSmartPointer.h>
 #include <vtkDoubleArray.h>
- 
- 
-void TriangulateTerrain();
+
+#include <vtkGraphicsFactory.h>
+#include <vtkImagingFactory.h>
+#include <vtkWindowToImageFilter.h>
+#include <vtkPNGWriter.h>
+
+typedef double (*Func3d)(double, double);
+
+void InitOffscreen(bool offscreen);
+vtkSmartPointer<vtkDelaunay2D> TriangulateTerrain(Func3d func);
+void PlotPoints(vtkSmartPointer<vtkDelaunay2D> delaunay, bool offscreen=true);
+
+double F(double x, double y)
+{
+  return atan(x*y*y);
+}
 
 int main (int argc, char ** argv)
 {
-	TriangulateTerrain();	
+	auto delaunay = TriangulateTerrain(&F);	
+  PlotPoints(delaunay,false);
 	return 0;
 }
-  
-void TriangulateTerrain()
+
+void InitOffscreen(bool offscreen)
 {
- 
+  // Graphics Factory
+  vtkSmartPointer<vtkGraphicsFactory> graphics_factory
+    = vtkGraphicsFactory::New();
+  graphics_factory->SetOffScreenOnlyMode( offscreen );
+  graphics_factory->SetUseMesaClasses( offscreen );
+
+  // Imaging Factory
+  vtkSmartPointer<vtkImagingFactory>  imaging_factory
+    = vtkImagingFactory::New();
+  imaging_factory->SetUseMesaClasses( offscreen );
+}
+
+vtkSmartPointer<vtkDelaunay2D> TriangulateTerrain(Func3d func)
+{
+
 	vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
   vtkSmartPointer<vtkDoubleArray> scalars = vtkSmartPointer<vtkDoubleArray>::New();
  
@@ -40,30 +68,37 @@
 	{
 		for(double y = -3; y < 3; y += 6.0/GridSize)
 		{
-      double z = sin(x*y);
+      double z = func(x,y);
 			vtkIdType idx = points->InsertNextPoint(x, y, z);
       scalars -> InsertTuple(idx, &z);
 		}
- 
 	}
  
 	//add the grid points to a polydata object
-	vtkPolyData* polydata = vtkPolyData::New();
+	vtkSmartPointer<vtkPolyData> polydata = vtkPolyData::New();
 	polydata->SetPoints(points);
   polydata->GetPointData() -> SetScalars(scalars);
  
 	//triangulate the grid points
-	vtkDelaunay2D* delaunay = vtkDelaunay2D::New();
+	vtkSmartPointer<vtkDelaunay2D> delaunay = vtkDelaunay2D::New();
 	delaunay->SetInput(polydata);
 	delaunay->Update();
 
+  return delaunay;
+}
+
+void PlotPoints(vtkSmartPointer<vtkDelaunay2D> delaunay, bool offscreen)
+{
+  InitOffscreen(offscreen);
+
+
   //Normals for Gourad shading
-  vtkPolyDataNormals* normals = vtkPolyDataNormals::New();
+  vtkSmartPointer<vtkPolyDataNormals> normals = vtkPolyDataNormals::New();
   normals -> SetInputConnection(delaunay ->GetOutputPort() );
   normals -> SetFeatureAngle(60);
 
   //Set the colours for the rendering
-  vtkLookupTable* lut = vtkLookupTable::New();
+  vtkSmartPointer<vtkLookupTable> lut = vtkLookupTable::New();
   lut -> SetHueRange(0.66667, 0.0);
   lut -> SetNumberOfColors(256);
   lut -> SetRampToLinear();
@@ -71,42 +106,65 @@
 
   
 	// map the contours to graphical primitives
-	vtkPolyDataMapper *contMapper = vtkPolyDataMapper::New();
+	vtkSmartPointer<vtkPolyDataMapper> contMapper = vtkPolyDataMapper::New();
   contMapper->SetInput(normals -> GetOutput() );
   contMapper->SetScalarRange(0,1);
   contMapper->SetLookupTable(lut);
   contMapper->ImmediateModeRenderingOn();
 
  
- 	 // create an actor for the contours
-	vtkActor *contActor = vtkActor::New();
+  // create an actor for the contours
+	vtkSmartPointer<vtkActor> contActor = vtkActor::New();
 	contActor->SetMapper(contMapper);
  
-	  // a renderer and render window
-	vtkRenderer *ren1 = vtkRenderer::New();
-	vtkRenderWindow *renWin = vtkRenderWindow::New();
+  // a renderer and render window
+	vtkSmartPointer<vtkRenderer> ren1 = vtkRenderer::New();
+	vtkSmartPointer<vtkRenderWindow> renWin = vtkRenderWindow::New();
+  if(offscreen)
+  {
+    renWin->SetOffScreenRendering(1);
+  }
 	renWin->AddRenderer(ren1);
- 
-  	// an interactor
-	vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
-	iren->SetRenderWindow(renWin);
-  
-  //The style
-  vtkInteractorStyleTerrain *terrain_style = vtkInteractorStyleTerrain::New();
- 
+
+  // an interactor
+  vtkSmartPointer<vtkRenderWindowInteractor> iren 
+    = vtkRenderWindowInteractor::New();
+  iren->SetRenderWindow(renWin);  
+
 	// add the actors to the scene
-	ren1->SetBackground(0.6,0.6,0.6); 
+	ren1->SetBackground(0.0,0.0,0.0);
   ren1->GetActiveCamera()->SetViewUp(0,0,1);
+  ren1->GetActiveCamera()->SetPosition(1, -1, -1);
   ren1->GetActiveCamera()->SetParallelProjection(1);
 	ren1->AddActor(contActor);
-  iren -> SetInteractorStyle(terrain_style);
   ren1->ResetCamera();
-
  
-  	// render an image (lights and cameras are created automatically)
+  // render an image (lights and cameras are created automatically)
 	renWin->Render();
+
+  if(offscreen)
+  {
+    // Write to file
+    vtkSmartPointer<vtkWindowToImageFilter> windowToImageFilter = 
+      vtkSmartPointer<vtkWindowToImageFilter>::New();
+    windowToImageFilter->SetInput(renWin);
+    windowToImageFilter->Update();
  
-  	// begin mouse interaction
-	iren->Start();
+    vtkSmartPointer<vtkPNGWriter> writer = vtkSmartPointer<vtkPNGWriter>::New();
+    writer->SetFileName("screenshot.png");
+    writer->SetInput(windowToImageFilter->GetOutput());
+    writer->Write();
+  }
+  else
+  {
+    //The style
+    vtkSmartPointer<vtkInteractorStyleTerrain> terrain_style 
+      = vtkInteractorStyleTerrain::New();
  
+    // Set the style
+    iren -> SetInteractorStyle(terrain_style);
+  
+    // begin mouse interaction
+    iren->Start(); 
+  }
 }