1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/experiments/torus/CMakeLists.txt	Wed Aug 04 00:10:25 2010 -0400
     1.3 @@ -0,0 +1,10 @@
     1.4 +set(QT_USE_QTOPENGL 1)
     1.5 +set(QT_USE_QTXML 1)
     1.6 +include(${QT_USE_FILE})
     1.7 +
     1.8 +add_executable(torus torus.cpp)
     1.9 +target_link_libraries(torus cloudy ${QT_LIBRARIES})
    1.10 +
    1.11 +add_executable(gentorus gentorus.cpp)
    1.12 +target_link_libraries(torus cloudy)
    1.13 +

     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/experiments/torus/gentorus.cpp	Wed Aug 04 00:10:25 2010 -0400
     2.3 @@ -0,0 +1,18 @@
     2.4 +#include "torus.hpp"
     2.5 +
     2.6 +int main()
     2.7 +{
     2.8 +  Torus torus(1.0, 0.2);
     2.9 +  Data_cloud cloud;
    2.10 +  size_t N = 10000;
    2.11 +  generate(torus, N, cloud);
    2.12 +
    2.13 +  for (size_t i = 0; i < cloud.size(); ++i)
    2.14 +    {
    2.15 +      const uvector &uv = cloud[i];
    2.16 +
    2.17 +      std::cout << uv(0) << " "
    2.18 +		<< uv(1) << " "
    2.19 +		<< uv(2) << "\n";
    2.20 +    }
    2.21 +}

     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/experiments/torus/pctviewer.hpp	Wed Aug 04 00:10:25 2010 -0400
     3.3 @@ -0,0 +1,54 @@
     3.4 +#include <cloudy/view/Viewer.hpp>
     3.5 +#include <cloudy/misc/Program_options.hpp>
     3.6 +
     3.7 +#include <qapplication.h>
     3.8 +#include <qmainwindow.h>
     3.9 +
    3.10 +#include <QDialog>
    3.11 +#include <QPushButton>
    3.12 +#include <QHBoxLayout>
    3.13 +
    3.14 +//using namespace cloudy::view;
    3.15 +
    3.16 +bool setup(cloudy::view::Viewer &w,
    3.17 +           const std::map<std::string, std::string> &options,
    3.18 +           const std::vector<std::string> &parameters);
    3.19 +
    3.20 +int main(int argc, char** argv)
    3.21 +{  
    3.22 +   std::map<std::string, std::string> options;
    3.23 +   std::vector<std::string> parameters;
    3.24 +   cloudy::misc::get_options(argc, argv, options, parameters);
    3.25 +
    3.26 +   QApplication app(argc, argv);
    3.27 +   QMainWindow main;
    3.28 +
    3.29 +   cloudy::view::Viewer w(0);
    3.30 +   
    3.31 +   if (!setup(w, options, parameters))
    3.32 +      return -1;
    3.33 +   
    3.34 +   w.setMinimumWidth(800);
    3.35 +   w.show();
    3.36 +
    3.37 +
    3.38 +   QDialog *dialog = new QDialog;
    3.39 +   QHBoxLayout *layout = new QHBoxLayout;
    3.40 +   QToolBox *tab = new QToolBox();
    3.41 +   layout->addWidget(tab);
    3.42 +   w.fill_dialog(tab);
    3.43 +   dialog->setMaximumWidth(250);
    3.44 +   dialog->setLayout(layout);
    3.45 +
    3.46 +   QWidget *centralWidget = new QWidget;
    3.47 +   QHBoxLayout *centralLayout = new QHBoxLayout;
    3.48 +   centralLayout->addWidget(&w);
    3.49 +   centralLayout->addWidget(dialog);
    3.50 +   centralWidget->setLayout(centralLayout);
    3.51 +
    3.52 +   main.setCentralWidget(centralWidget);
    3.53 +   main.resize(1050, 600);
    3.54 +   main.show(); 
    3.55 +
    3.56 +   return app.exec();
    3.57 +}

     4.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     4.2 +++ b/experiments/torus/torus.cpp	Wed Aug 04 00:10:25 2010 -0400
     4.3 @@ -0,0 +1,85 @@
     4.4 +#include <cloudy/linear/Covariance.hpp>
     4.5 +#include <cloudy/misc/Progress.hpp>
     4.6 +
     4.7 +#include "torus.hpp"
     4.8 +#include "pctviewer.hpp"
     4.9 +
    4.10 +using namespace cloudy::view;
    4.11 +using namespace cloudy;
    4.12 +
    4.13 +void diagonalize(const Data_cloud_ptr &covariance,
    4.14 +		 Data_cloud_ptr &normals,
    4.15 +		 Data_cloud_ptr &K1,
    4.16 +		 Data_cloud_ptr &K2,
    4.17 +		 Scalar_field_ptr &anisotropy)
    4.18 +{
    4.19 +  cloudy::misc::Progress_display progress(covariance->size(), std::cerr);
    4.20 +   for (Data_cloud::iterator it = covariance->begin();
    4.21 +        it != covariance->end(); ++it)
    4.22 +   {
    4.23 +      std::vector<double> V;
    4.24 +      std::vector<uvector> D;
    4.25 +
    4.26 +      linear::covariance_extract_eigen
    4.27 +	(linear::matrix_from_covariance_3(*it), V, D);
    4.28 +      linear::covariance_sort_eigen(V, D);
    4.29 +      normals->push_back(D[0]);
    4.30 +      K1->push_back(D[1]);
    4.31 +      K2->push_back(D[2]);
    4.32 +      anisotropy->push_back(V[1]/(V[0] + V[1] + V[2]));
    4.33 +      ++progress;
    4.34 +   }
    4.35 +   std::cerr << "done\n";
    4.36 +}
    4.37 +
    4.38 +
    4.39 +bool setup(cloudy::view::Viewer &w,
    4.40 +           const std::map<std::string, std::string> &options,
    4.41 +           const std::vector<std::string> &parameters)
    4.42 +{
    4.43 +  using namespace cloudy::view;
    4.44 +
    4.45 +  Torus torus(1.0, 0.2);
    4.46 +
    4.47 +  Data_cloud_ptr cloud (new Data_cloud());
    4.48 +  Data_cloud_ptr covariance (new Data_cloud());
    4.49 +
    4.50 +  assert(parameters.size() == 2);
    4.51 +
    4.52 +  std::ifstream is(parameters[0].c_str());
    4.53 +  std::ifstream isP(parameters[1].c_str());
    4.54 +  
    4.55 +  std::cerr << "loading " << parameters[0] << ".. ";
    4.56 +  cloudy::load_cloud(is, *cloud);
    4.57 +  std::cerr << "done\n";
    4.58 +  
    4.59 +  std::cerr << "loading " << parameters[1] << ".. ";
    4.60 +  cloudy::load_cloud(isP, *covariance);
    4.61 +
    4.62 +  assert(cloud->size() == covariance->size());
    4.63 +
    4.64 +  std::cerr << "diagonalizing covariance matrices..\n";
    4.65 +  Data_cloud_ptr normals (new Data_cloud());
    4.66 +  Data_cloud_ptr K1 (new Data_cloud());
    4.67 +  Data_cloud_ptr K2 (new Data_cloud());
    4.68 +  Scalar_field_ptr anisotropy (new Scalar_field());
    4.69 +  diagonalize(covariance, normals, K1, K2, anisotropy);
    4.70 +
    4.71 +  Data_cloud_ptr rnormals (new Data_cloud());
    4.72 +  double total_alpha (0.0);
    4.73 +  for (size_t i = 0; i < cloud->size(); ++i)
    4.74 +    {
    4.75 +      rnormals->push_back(torus.get_normal((*cloud)[i]));
    4.76 +      double alpha = acos( fabs(ublas::inner_prod((*rnormals)[i], (*normals)[i])) );
    4.77 +      total_alpha += alpha;
    4.78 +      //      std::cerr << alpha << "\n";
    4.79 +    }  
    4.80 +  std::cerr << "Average error: " << total_alpha/(cloud->size()) << "\n";
    4.81 +  
    4.82 +  Cloud_drawer *c = new Cloud_drawer("Cloud", cloud);
    4.83 +  w.add_drawer(Drawer_ptr(c));
    4.84 +  w.add_drawer(Drawer_ptr(new Direction_drawer("Normals", *c, normals)));
    4.85 +  w.add_drawer(Drawer_ptr(new Direction_drawer("Real Normals", *c, rnormals)));
    4.86 +   
    4.87 +   return true;
    4.88 +}

     5.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.2 +++ b/experiments/torus/torus.hpp	Wed Aug 04 00:10:25 2010 -0400
     5.3 @@ -0,0 +1,110 @@
     5.4 +#ifndef TORUS_HPP 
     5.5 +#define TORUS_HPP
     5.6 +
     5.7 +#include <cloudy/random/Random.hpp>
     5.8 +#include <fstream>
     5.9 +#include <cloudy/Cloud.hpp>
    5.10 +
    5.11 +using namespace cloudy;
    5.12 +inline uvector uv3 (double x, double y, double z)
    5.13 +{
    5.14 +  uvector v(3);
    5.15 +  v(0) = x;
    5.16 +  v(1) = y;
    5.17 +  v(2) = z;
    5.18 +  return v;
    5.19 +}
    5.20 +
    5.21 +inline uvector rotate_around_z(const uvector &v, double theta)
    5.22 +{
    5.23 +  return uv3(cos(theta)*v(0) - sin(theta)*v(1),
    5.24 +	     sin(theta)*v(0) + cos(theta)*v(1),
    5.25 +	     v(2));
    5.26 +}
    5.27 +
    5.28 +enum
    5.29 +  {
    5.30 +    // Angle in (x,y), corresponding to radius R: theta
    5.31 +    COORD_THETA,
    5.32 +    COORD_PHI,
    5.33 +    COORD_L
    5.34 +  };
    5.35 +
    5.36 +#define DBG(x) std::cerr << "in " <<  __FUNCTION__ << ", "<< #x << " = " << x << "\n"
    5.37 +
    5.38 +struct Torus
    5.39 +{
    5.40 +  double _r,  _R; 
    5.41 +
    5.42 +  uvector xyz_to_torus (const uvector &uv) const
    5.43 +  {
    5.44 +    double theta = atan2(uv(1), uv(0));
    5.45 +    uvector uv_in_xz = rotate_around_z(uv, -theta);
    5.46 +    uv_in_xz(0) -= _R;
    5.47 +    // uv_in_xz lies in the plan (x,z).
    5.48 +
    5.49 +    double phi = atan2(uv_in_xz(2), uv_in_xz(0));
    5.50 +    double l = ublas::norm_2(uv_in_xz);
    5.51 +
    5.52 +    //DBG(uv); DBG(uv_in_xz);
    5.53 +    // DBG(theta);  DBG(phi);  DBG(l); 
    5.54 +
    5.55 +    return uv3(theta, phi, l);
    5.56 +  }
    5.57 +
    5.58 +  uvector
    5.59 +  torus_to_xyz (const uvector &tor) const
    5.60 +  {
    5.61 +    double theta = tor(COORD_THETA), phi = tor(COORD_PHI),
    5.62 +      L = tor(COORD_L);
    5.63 +    
    5.64 +    uvector uv_in_xz = uv3(L * cos(phi) + _R, 0, L * sin(phi));
    5.65 +    return rotate_around_z(uv_in_xz, theta);
    5.66 +  }
    5.67 +
    5.68 +  uvector
    5.69 +  torus_to_normal (const uvector &tor) const
    5.70 +  {
    5.71 +    double theta = tor(COORD_THETA), phi = tor(COORD_PHI);
    5.72 +    
    5.73 +    uvector uv_in_xz = uv3(cos(phi), 0, sin(phi));
    5.74 +    return rotate_around_z(uv_in_xz, theta);
    5.75 +  }
    5.76 +
    5.77 +public:
    5.78 +  Torus(double R, double r): _r(r), _R(R) {}
    5.79 +
    5.80 +  bool is_inside(const uvector &uv) const
    5.81 +  {
    5.82 +    uvector tor = xyz_to_torus(uv);
    5.83 +    return (tor(COORD_L) <= _r);
    5.84 +  }
    5.85 +
    5.86 +  uvector project(const uvector &uv) const
    5.87 +  {
    5.88 +    uvector tor = xyz_to_torus(uv);
    5.89 +    tor(COORD_L) = _r;
    5.90 +    return torus_to_xyz(tor);
    5.91 +  }
    5.92 +
    5.93 +  uvector get_normal(const uvector &uv) const  
    5.94 +  {
    5.95 +    return torus_to_normal(xyz_to_torus(uv));
    5.96 +  }
    5.97 +};
    5.98 +
    5.99 +inline void generate(const Torus &torus, size_t N,
   5.100 +		     Data_cloud &cloud)
   5.101 +{
   5.102 +  random::Random_vector_in_ball<uvector> randball (3, torus._R + torus._r);
   5.103 +  boost::mt19937 engine;
   5.104 +
   5.105 +  while (cloud.size() != N)
   5.106 +    {
   5.107 +      uvector uv = randball(engine);
   5.108 +      if (torus.is_inside(uv))
   5.109 +	cloud.push_back(torus.project(uv));
   5.110 +    }
   5.111 +}
   5.112 +
   5.113 +#endif