Mercurial > hg > kwantix
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(); + } }