Commits

Anonymous committed e22ca29

1. add reassigning ids in ContainerOfIdItems. 2. add generating uniform poits on SurfaceMesh

  • Participants
  • Parent commits 812dbd9

Comments (0)

Files changed (13)

+syntax: regexp
+(^|/)\.DS_Store$
+^\.idea$
+\.sublime-
+^nbproject$
+^build$
+^docs/_build$
+^dist$
+^[^/]*.egg-info$
+^\.settings$
+~$
+\.bak$
+\.orig$
+\.pyc$
+\.pyo$
+\.swp$
+\.tmp$
+^_generated_
+^\.cache_
+desktop\.ini$

cmake/FindSlog.cmake

+# - Find Slog (C++ logging)
+# This module defines
+# Slog_INCLUDE_DIR, where to find Slog headers
+# Slog_LIB, Slog libraries
+# Slog_FOUND, If false, do not try to use Slog
+
+# Consider as Idontcare licence.
+
+find_path(Slog_INCLUDE_DIR Slog/Logger.h PATHS
+            /usr/local/include
+            /opt/local/include)
+find_library(Slog_LIB NAMES slog PATHS
+            /usr/local/lib
+            /opt/local/lib)
+
+if (Slog_LIB AND Slog_INCLUDE_DIR)
+    set(Slog_FOUND TRUE)
+    set(Slog_LIBS ${Slog_LIB})
+    message(STATUS "found Slog: ${Slog_LIB}")
+else ()
+    set(Slog_FOUND FALSE)
+    if (Slog_FIND_REQUIRED)
+        message(FATAL_ERROR "Slog NOT found")
+    else ()
+        message(STATUS "Slog not found")
+    endif ()
+endif ()
+
+mark_as_advanced(
+    Slog_LIB
+    Slog_INCLUDE_DIR
+  )
     class BaseId
     {
         public:
-            BaseId(){eventCounter = 0;};
+            BaseId(){id = 0; eventCounter = 0;};
             ~BaseId(){};
             void setId(size_t id) {this->id = id;};
             size_t getId() const {return id;};

include/ContainerOfIdItems.h

 #include <string>
 #include <memory>
 #include <stdexcept>
+#include <Slog/Logger.h>
 #ifdef HAVE_KLUBOK
 #include <Klubok/thread.h>
 #endif
     using std::string;
     using std::auto_ptr;
     using std::invalid_argument;
+    using slog::Logger;
 #ifdef HAVE_KLUBOK
     using klubok::Mutex;
 #endif
             void add(vector<ItemType* > &items);
             void add(const ContainerOfIdItems<ItemType> &other);
             template<class ContainerType>
-            void add(const vector<const ContainerType*> &others);
+            void add(const vector<ContainerType*> &others);
             void transfer(ContainerOfIdItems<ItemType> &other);
             template<class ContainerType>
             void transfer(vector<ContainerType*> &others);
             
+            void assignIncrementalIds();
             auto_ptr<ItemType> remove(size_t index);
             void drop();
             void clear();
     // =========================================================================
     template<class ItemType>
     template<class ContainerType>
-    void ContainerOfIdItems<ItemType>::add(const vector<const ContainerType*>
+    void ContainerOfIdItems<ItemType>::add(const vector<ContainerType*>
                                                                         &others)
     {
         for (size_t i = 0; i < others.size(); i++)
     template<class ItemType>
     void ContainerOfIdItems<ItemType>::getIds(vector<size_t> &res) const
     {
+        Logger::getLogger().getDebugSink()<<"ContainerOfIdItems make vector "
+                                          <<"of all ids, size "<<items.size();
         res.resize(0);
         res.reserve(items.size());
         for (size_t i = 0; i < items.size(); i++)
-        {
             res.push_back(items[i]->getId());
-        }
     }
 
     // =========================================================================
                                              const vector<size_t> &ids,
                                              vector<size_t> &indices) const
     {
+        Logger &logger = Logger::getLogger();
         if (ids.empty())
+        {
+            logger.getDebugSink()<<" no ids provided, no processing required";
             return;
+        }
         //create back map
+        logger.getDebugSink()<<"get all indices";
         vector<size_t> allids, backMap;
         getIds(allids);
+        logger.getDebugSink()<<"create back map id -> index";
+        logger.getDebugSink()<<"find maximal id";
         size_t max = allids[0];
         for (size_t i = 1; i < allids.size(); i++)
         {
             if (allids[i] > max)
                 max = allids[i];
         }
+        logger.getDebugSink()<<"maximal id = "<<max
+                             <<" resize back map vector to max+1";
         backMap.resize(max+1, items.size());
         for (size_t i = 0; i < allids.size(); i++)
             backMap[allids[i]] = i;
 
         //process requested ids
+        logger.getDebugSink()<<"iterate over ids and set corresponding indices, "
+                             <<" size "<<ids.size();
         indices.resize(ids.size(), 0);
         for (size_t i = 0; i < ids.size(); i++)
         {
     void ContainerOfIdItems<ItemType>::getByIds(const vector<size_t> &ids,
                                          vector<ItemType*> &result) const
     {
+        Logger::getLogger().getDebugSink()<<"ContainerOfIdItems get items "
+                                    <<"specified by their ids "<<std::hex<<this;
         vector<size_t> indices;
+        Logger::getLogger().getDebugSink()<<"find indices of items";
         getIndicesByIds(ids, indices);
+        Logger::getLogger().getDebugSink()<<"resize and fill vector of items, "
+                                          <<"size "<<indices.size();
         result.resize(indices.size());
         for (size_t i = 0; i < indices.size(); i++)
-        {
             result[i] = items[indices[i]];
-        }
     }
     // =========================================================================
     /// @brief Return list of items with given ids.
         for (size_t i = 0; i < ids.size(); i++)
             backMap[ids[i]] = i;
     }
+
+    // =========================================================================
+    /// @brief Assigns incremental ids to all items in the container.
+    //
+    // Old ids are dropped and items are assigned an id according to
+    // their position in the container.
+    // =========================================================================
+    template<class ItemType>
+    void ContainerOfIdItems<ItemType>::assignIncrementalIds()
+    {
+        for (size_t i = 0; i < items.size(); i++)
+            items[i]->setId(i);
+    }
 } //namespace meshlib 
 #endif /* end of include guard: CONTAINEROFIDITEMS_H_JOPBB5RI */

include/RandomGenerator.h

 #include <gsl/gsl_rng.h>
 #endif     // -----  HAVE_GSL  -----
 #include <ostream>
+#include <vector>
 
 namespace meshlib
 {
+    using std::vector;
     
     // =========================================================================
     /// @brief Wrapper around random generators.
             ~RandomGenerator();
             void init_rand();
             double get_rand() const;
+            void sampleBySize(size_t numSamples,
+                              const vector<double> &sizes,
+                              vector<size_t> &inds) const;
             friend std::ostream &operator<<(std::ostream &stream,
                                                 const RandomGenerator &object);
         private:

include/SurfaceMesh.h

             double getArea() const;
             double getArea(vector<double> &res) const;
             bool containsPolygon(const AbstractPolygon *polygon) const;
+            void getUniformRandomPoints(const RandomGenerator &generator,
+                                        int numPoints, 
+                                        vector<Point3D > &res) const;
             void getRandomPointsFixN(const RandomGenerator &generator,
                     int numPoints, vector<vector<Point3D> > &res) const;
             pair<double, const AbstractPolygon*>

service/MeshSrvBase.h

 #ifndef MESHSRVBASE_JP3C62R
 #define MESHSRVBASE_JP3C62R
 
+#include "ContainerOfIdItems.h"
 #include <LinAlg/Vector3D.h>
 #include <Slog/Logger.h>
 #include <vector>
 namespace meshlib
 { 
     class SurfaceMesh;
+    class SurfaceMeshSet;
+    class VolumeMesh;
 } //namespace meshlib
 
 namespace meshsrv
     using slog::Logger;
     using std::vector;
     using std::string;
+    using meshlib::ContainerOfIdItems;
     using meshlib::SurfaceMesh;
+    using meshlib::SurfaceMeshSet;
+    using meshlib::VolumeMesh;
     using linalg::Vector3D;
     
     // =========================================================================
         public:
             MeshSrvBase() : logger(Logger::getLogger()) {};
             ~MeshSrvBase(){};
-
-            void getPolygonNormals(const vector<int> &polygonIds,
+            
+            //void dropDegeneratedPolygons(const SurfaceMesh &inMesh,
+                                         //SurfaceMesh &outMesh);
+            void refineMesh(const SurfaceMesh &inMesh,
+                            double threshold,
+                            SurfaceMesh &outMesh);
+            void triangulateMesh(const SurfaceMesh &inMesh,
+                                 SurfaceMesh &outMesh);
+            void combineMeshes(const SurfaceMeshSet &meshSet,
+                               SurfaceMesh &mesh);
+            void getPolygonNormals(const vector<size_t> &polygonIds,
                                    const SurfaceMesh &mesh,
                                    vector<Vector3D> &normals);
-            void getPolygonCenters(const vector<int> &polygonIds,
+            void getPolygonCenters(const vector<size_t> &polygonIds,
                                    const SurfaceMesh &mesh,
+                                   double relTol, double absTol,
                                    vector<Vector3D> &centers);
-            void getPolygonAreas(const vector<int> &polygonIds,
+            void getPolygonAreas(const vector<size_t> &polygonIds,
                                  const SurfaceMesh &mesh,
                                  vector<double> &areas);
             void getMeshArea(const SurfaceMesh &mesh,
                              double &area);
+
+            void getPolyhedronCenters(const vector<size_t> &polyhedronIds,
+                                      const VolumeMesh &mesh,
+                                      vector<Vector3D> &centers);
+            void getPolyhedronSurfaceAreas(const vector<size_t> &polyhedronIds,
+                                           const VolumeMesh &mesh,
+                                           vector<double> &areas);
+            void getPolyhedronVolumes(const vector<size_t> &polyhedronIds,
+                                      const VolumeMesh &mesh,
+                                      vector<double> &volumes);
+            void getMeshVolume(const VolumeMesh &mesh,
+                               double &volume);
+
+           //void getUniformRandomPoints(size_t numPoints,
+                                       //const SurfaceMesh &mesh,
+                                       //::geometry::Points3D& _return,
+           //void getRandomPointsForPolygon(std::vector< ::geometry::Points3D> & _return, const int32_t numPoints, const std::vector<int32_t> & polygonIds, const SurfaceMeshWrap& mesh) {
+
+            void getBuildInfo(string& info);
+
         private:
             Logger &logger;
-            bool checkMeshEmpty(const SurfaceMesh &mesh) const;
-            bool checkVectorEmpty(const vector<int> &vec,
+            template <class T>
+            bool checkContainerEmpty(const ContainerOfIdItems<T> &mesh,
+                                     const string &contName) const;
+            template <class T>
+            bool checkVectorEmpty(const vector<T> &vec,
                                   const string &vecName) const;
     };
+    //
+    // =========================================================================
+    /// @brief Check if vector is empty, use vector name to log the result.
+    // =========================================================================
+    template <class T>
+    bool MeshSrvBase::checkVectorEmpty(const vector<T> &vec,
+                                       const string &vecName) const
+    {
+        if (vec.empty())
+        {
+            logger.getDebugSink()<<"empty vector "<<vecName;
+            return true;
+        }
+        return false;
+    }
 
+    template <class T>
+    bool MeshSrvBase::checkContainerEmpty(const ContainerOfIdItems<T> &cont,
+                                          const string &contName) const
+    {
+        if (cont.getNumItems() == 0)
+        {
+            logger.getDebugSink()<<"empty "<<contName;
+            return true;
+        }
+        return false;
+    }
 
-//void getVoxelSurfaceAreas(std::vector<double> & _return, const std::vector<int32_t> & voxelIds, const VolumeMeshWrap& mesh) {
-  //// Your implementation goes here
-  //printf("getVoxelSurfaceAreas\n");
-//}
 
-//void getVoxelVolumes(std::vector<double> & _return, const std::vector<int32_t> & voxelIds, const VolumeMeshWrap& mesh) {
-  //// Your implementation goes here
-  //printf("getVoxelVolumes\n");
-//}
-
-//void getVoxelCenters(std::vector< ::geometry::Vector> & _return, const std::vector<int32_t> & voxelIds, const VolumeMeshWrap& mesh) {
-  //// Your implementation goes here
-  //printf("getVoxelCenters\n");
-//}
-
-//double getMeshVolume(const VolumeMeshWrap& mesh) {
-  //// Your implementation goes here
-  //printf("getMeshVolume\n");
-//}
 
 
     ////  void intersectSegmentMesh(HitPointsAndPolyLoads& _return, const  ::geometry::Points3D& starts, const  ::geometry::Points3D& stops, const MeshSet& meshSet) {
     ////    printf("intersectMeshPhiPlane\n");
     ////  }
     ////
-    ////  void combineMeshes( ::componentsdb2::SurfaceMesh& _return, const MeshSet& meshSet) {
-    ////    // Your implementation goes here
-    ////    printf("combineMeshes\n");
-    ////  }
-    ////
-    ////  void triangulateMesh( ::componentsdb2::SurfaceMesh& _return, const SurfaceMeshWrap& mesh) {
-    ////    // Your implementation goes here
-    ////    printf("triangulateMesh\n");
-    ////  }
-    ////
-    ////  void refineMesh( ::componentsdb2::SurfaceMesh& _return, const SurfaceMeshWrap& mesh, const double threshold) {
-    ////    // Your implementation goes here
-    ////    printf("refineMesh\n");
-    ////  }
-    ////
-    ////  void dropDegeneratedPolygons( ::componentsdb2::SurfaceMesh& _return, const SurfaceMeshWrap& mesh) {
-    ////    // Your implementation goes here
-    ////    printf("dropDegeneratedPolygons\n");
-    ////  }
     ////
     ////  void transformPoints( ::geometry::Points3D& _return, const  ::geometry::Points3D& points, const  ::geometry::Transformation& trans) {
     ////    // Your implementation goes here
     ////  }
     ////
     ////
-    ////  void getPolygonCenters( ::geometry::Points3D& _return, const std::vector<int32_t> & polygonIds, const SurfaceMeshWrap& mesh) {
-    ////    // Your implementation goes here
-    ////    printf("getPolygonCenters\n");
-    ////  }
-    ////
-    ////  void getPolygonNormals(std::vector< ::geometry::Vector> & _return, const std::vector<int32_t> & polygonIds, const SurfaceMeshWrap& mesh) {
-    ////    // Your implementation goes here
-    ////    printf("getPolygonNormals\n");
-    ////  }
-    ////
-    ////  double getMeshArea(const SurfaceMeshWrap& mesh) {
-    ////    // Your implementation goes here
-    ////    printf("getMeshArea\n");
-    ////  }
-    ////
-    ////  void getVoxelSurfaceAreas(std::vector<double> & _return, const std::vector<int32_t> & voxelIds, const VolumeMeshWrap& mesh) {
-    ////    // Your implementation goes here
-    ////    printf("getVoxelSurfaceAreas\n");
-    ////  }
-    ////
-    ////  void getVoxelVolumes(std::vector<double> & _return, const std::vector<int32_t> & voxelIds, const VolumeMeshWrap& mesh) {
-    ////    // Your implementation goes here
-    ////    printf("getVoxelVolumes\n");
-    ////  }
-    ////
-    ////  void getVoxelCenters(std::vector< ::geometry::Vector> & _return, const std::vector<int32_t> & voxelIds, const VolumeMeshWrap& mesh) {
-    ////    // Your implementation goes here
-    ////    printf("getVoxelCenters\n");
-    ////  }
-    ////
-    ////  double getMeshVolume(const VolumeMeshWrap& mesh) {
-    ////    // Your implementation goes here
-    ////    printf("getMeshVolume\n");
-    ////  }
-    ////
-    ////  void getUniformRandomPoints( ::geometry::Points3D& _return, const int32_t numPoints, const SurfaceMeshWrap& mesh) {
-    ////    // Your implementation goes here
-    ////    printf("getUniformRandomPoints\n");
-    ////  }
-    ////
-    ////  void getRandomPointsForPolygon(std::vector< ::geometry::Points3D> & _return, const int32_t numPoints, const std::vector<int32_t> & polygonIds, const SurfaceMeshWrap& mesh) {
-    ////    // Your implementation goes here
-    ////    printf("getRandomPointsForPolygon\n");
-    ////  }
-    ////
-    ////  void getBuildInfo(std::string& _return) {
-    ////    // Your implementation goes here
-    ////    printf("getBuildInfo\n");
-    ////  }
     ////
     ////  int64_t getTimeStamp() {
     ////    // Your implementation goes here

src/RandomGenerator.cpp

 // Licensed under GPLv3 or later, see the COPYING file.
 
 #include "RandomGenerator.h"
+#include <Slog/Logger.h>
+#include <stdexcept>
 #include <ctime>
+#include <algorithm>
 
+using namespace std;
 using namespace meshlib;  
+using slog::Logger;
 // =========================================================================
 /// @brief Constructor
 ///
 #endif     // -----  not HAVE_GSL  -----
     return stream;
 }
+
+// =========================================================================
+/// @brief Do uniform sampling of intervals according to their sizes.
+/// 
+//  Every interval is described by its size. The input size vector is not
+//  sorted. This method samples intervals uniformly, i.e. the probability
+//  of hitting the interval is proportional to its size. This is an
+//  example of reverse sampling theorem.
+//
+/// @param numSamples Number of samples to generate.
+/// @param sizes Vector of interval sizes.
+/// @param inds Indices of chosen intervals.
+// =========================================================================
+void RandomGenerator::sampleBySize(size_t numSamples,
+                                   const vector<double> &sizes,
+                                   vector<size_t> &inds) const
+{
+    Logger &logger = Logger::getLogger();
+    logger.getDebugSink()<<"RandomGenerator sample intervals by their size, "
+                         <<numSamples<<" samples to generate";
+    if (numSamples == 0)
+    {
+        logger.getDebugSink()<<"0 samples required, return";
+        inds.resize(0);
+        return;
+    }
+    logger.getDebugSink()<<"generate cumulative sum vector";
+    vector<double> cumsum(sizes.size(), 0);
+    cumsum[0] = sizes[0];
+    for (size_t i = 1; i < sizes.size(); i++)
+    {
+        if (sizes[i] == 0)
+        {
+            logger.getWarnSink()<<"interval "<<i
+                <<"has zero size, it will get no samples";
+        }
+        cumsum[i] = cumsum[i-1] + sizes[i];
+    }
+    double total = cumsum[sizes.size()-1];
+    logger.getDebugSink()<<"total sum "<<total;
+    if (total == 0)
+    {
+        logger.getErrorSink()<<"total sum is 0, no sampling will be done";
+        throw invalid_argument("Zero total sum of all interval sizes, "
+                               "no sampling could bed done");
+    }
+    logger.getDebugSink()<<"normalize cumulative sums";
+    for (size_t i = 0; i < cumsum.size(); i++)
+        cumsum[i] /= total;
+    logger.getDebugSink()<<"draw samples";
+    inds.resize(numSamples);
+    for (size_t i = 0; i < numSamples; i++)
+    {
+        double v = get_rand();
+        //see where the value fits into
+        vector<double>::const_iterator it = lower_bound(cumsum.begin(),
+                                                        cumsum.end(), v);
+        inds[i] = it - cumsum.begin();
+    }
+}

src/SurfaceMesh.cpp

 #include "BoundingBox.h"
 #include "Triangle.h"
 #include "Quadrangle.h"
+#include <Slog/Logger.h>
 #include <stdexcept>
 #include <sstream>
 
 using std::make_pair;
 using std::getline;
 using std::istringstream;
+using slog::Logger;
     
 // =========================================================================
 /// @brief Construct from a set of polygons.
 }
 
 // =========================================================================
+/// @brief Get random points uniformly on the surface.
+/// 
+//
+// Generate points distributed uniformly on the surface, i.e.
+// choose polygons according to their area/surface area.
+// =========================================================================
+void SurfaceMesh::getUniformRandomPoints(const RandomGenerator &generator,
+                            int numPoints, vector<Point3D > &res) const
+{
+    Logger &logger = Logger::getLogger();
+    logger.getDebugSink()<<"SurfaceMesh generate uniform random points"
+                         <<std::hex<<this;
+    logger.getDebugSink()<<"make vector of areas of all polygons, size "
+                         <<getNumItems();
+    vector<double> areas;
+    getArea(areas);
+    logger.getDebugSink()<<"choose polygons by their area, "
+                         <<numPoints<<" indices to get";
+    vector<size_t> inds;
+    generator.sampleBySize(numPoints, areas, inds);
+    logger.getDebugSink()<<"generate one point per chose polygon";
+    res.resize(inds.size());
+    for (size_t i = 0; i < inds.size(); i++)
+        items[inds[i]]->getRandomPoint(generator, res[i]);
+}
+
+// =========================================================================
 /// @brief Get minimal distance from point to component.
 /// 
 /// @param point Point to get distance for.

tests/container.cpp

             BOOST_REQUIRE_EQUAL(bm[i], 3);
     }
 }
+
+BOOST_AUTO_TEST_CASE(assignIncrementalIds)
+{
+    vector<size_t> ids;
+    cont1.getIds(ids);
+    BOOST_REQUIRE_EQUAL(ids[0], 5);
+    BOOST_REQUIRE_EQUAL(ids[1], 6);
+    BOOST_REQUIRE_EQUAL(ids[2], 4);
+    BOOST_REQUIRE_EQUAL(ids[3], 0);
+    BOOST_REQUIRE_EQUAL(ids[4], 2);
+    cont1.assignIncrementalIds();
+    cont1.getIds(ids);
+    BOOST_REQUIRE_EQUAL(ids[0], 0);
+    BOOST_REQUIRE_EQUAL(ids[1], 1);
+    BOOST_REQUIRE_EQUAL(ids[2], 2);
+    BOOST_REQUIRE_EQUAL(ids[3], 3);
+    BOOST_REQUIRE_EQUAL(ids[4], 4);
+}
 BOOST_AUTO_TEST_SUITE_END()
         MyConfig()
         {
             Logger &logger = Logger::getLogger();
-            //logger.setDebugSink(auto_ptr<Sink>(new CoutSink()));
+            logger.setDebugSink(auto_ptr<Sink>(new CoutSink()));
             logger.setInfoSink(auto_ptr<Sink>(new CoutSink()));
             logger.setWarnSink(auto_ptr<Sink>(new CoutSink()));
             logger.setErrorSink(auto_ptr<Sink>(new CoutSink()));
 #define BOOST_TEST_DYN_LINK
 #include <boost/test/unit_test.hpp>
 #include "RandomGenerator.h"
+#include <Slog/Logger.h>
+#include <cmath>
 #include <iostream>
 
 using namespace std;
 using namespace meshlib;
+using namespace slog;
 
 //tests for RandomGenerator
 BOOST_AUTO_TEST_SUITE(Random_Tests)
         BOOST_REQUIRE( (x<1) && (x>=0) );
     }
 }
+
+BOOST_AUTO_TEST_CASE(sampleBySize)
+{
+    Logger &logger = Logger::getLogger();
+    RandomGenerator gen;
+    vector<double> sizes(4, 0);
+    sizes[0] = .1*3;
+    sizes[1] = .1*3;
+    sizes[2] = .2*3;
+    sizes[3] = .6*3;
+    vector<size_t> inds;
+    gen.sampleBySize(10000, sizes, inds);
+
+    vector<size_t> freq(4, 0);
+    for (size_t i = 0; i < inds.size(); i++)
+        freq[inds[i]]++;
+    for (size_t i = 0; i < freq.size(); i++)
+    {
+        double generated = double(freq[i])/inds.size();
+        double ideal = double(sizes[i])/3.0;
+        BOOST_REQUIRE(fabs(generated-ideal)<.01);
+    }
+}
+
 BOOST_AUTO_TEST_SUITE_END()

tests/surfacemesh.cpp

 #include <boost/test/floating_point_comparison.hpp>
 #include "SurfaceMesh.h"
 #include "Triangle.h"
+#include "Quadrangle.h"
 #include <limits>
 #include <iostream>
 using namespace meshlib;
     BOOST_REQUIRE_EQUAL(m2.get().size(), 3);
     BOOST_REQUIRE_EQUAL(m2.getArea(), m.getArea());
 }
+
+BOOST_AUTO_TEST_CASE(getUniformRandomPoints)
+{
+    Point3D p1, p2, p3, p4;
+    p1[0] = -1; p1[1] = -1;
+    p2[0] = 1; p2[1] = -1;
+    p3[0] = 1; p3[1] = 1;
+    p4[0] = -1; p4[1] = 1;
+    Quadrangle q1(p1,p2,p3,p4);
+    p1[1] = 1; p2[1] = 1; p3[1] = 2; p4[1] = 2;
+    Quadrangle q2(p1,p2,p3,p4);
+    p1[1] = 2; p2[1] = 2; p3[1] = 2.5; p4[1] = 2.5;
+    Quadrangle q3(p1,p2,p3,p4);
+    SurfaceMesh mesh;
+    mesh.add(q1);
+    mesh.add(q2);
+    mesh.add(q3);
+    vector<Point3D> points;
+    RandomGenerator generator;
+    mesh.getUniformRandomPoints(generator, 10000,  points);
+    BOOST_REQUIRE_EQUAL(points.size(), 10000);
+    double f1 = 0, f2 = 0, f3 = 0;
+    for (size_t i = 0; i < points.size(); i++)
+    {
+        if (points[i][1] <1)
+            f1++;
+        else if (points[i][1]<2)
+            f2++;
+        else
+            f3++;
+    }
+    f1 /= points.size(); 
+    f2 /= points.size(); 
+    f3 /= points.size(); 
+    BOOST_REQUIRE(fabs(f1-4/7.)<.008);
+    BOOST_REQUIRE(fabs(f2-2/7.)<.008);
+    BOOST_REQUIRE(fabs(f3-1/7.)<.008);
+}
 BOOST_AUTO_TEST_SUITE_END()