Index: src/Shapes/BaseShapes.cpp
===================================================================
--- src/Shapes/BaseShapes.cpp	(revision 0d02fbf6b0433481fd17dc075f5b57ba0c694fe6)
+++ src/Shapes/BaseShapes.cpp	(revision c5186e909d411ff3c5b6bc38edd607738818eab0)
@@ -6,6 +6,9 @@
  */
 
+#include "Helpers/Assert.hpp"
 #include "Shapes/BaseShapes.hpp"
 #include "Shapes/BaseShapes_impl.hpp"
+
+#include <cmath>
 
 #include "LinearAlgebra/Vector.hpp"
@@ -14,4 +17,33 @@
   return point.NormSquared()<=1;
 }
+
+
+/**
+ * algorithm taken from http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere
+ * \param N number of points on surface
+ */
+std::vector<Vector> Sphere_impl::getHomogeneousPointsOnSurface(const int N) const {
+  std::vector<Vector> PointsOnSurface;
+
+  const double dlength = M_PI*(3.-sqrt(5.));
+  double length = 0;
+  const double dz = 2.0/N;
+  double z = 1. - dz/2.;
+  Vector point;
+  for (int ka = 0; ka<N; ka++){
+    const double r = sqrt(1.-z*z);
+    point.Zero();
+    point[0] = cos(length)*r;
+    point[1] = sin(length)*r;
+    point[2] = z;
+    PointsOnSurface.push_back(point);
+    z = z - dz;
+    length = length + dlength;
+  }
+
+  ASSERT(PointsOnSurface.size() == N, "Sphere_impl::getHomogeneousPointsOnSurface() did not create enough points.");
+  return PointsOnSurface;
+}
+
 
 Shape Sphere(){
@@ -24,4 +56,13 @@
 }
 
+/**
+ * \param N number of points on surface
+ */
+std::vector<Vector> Cuboid_impl::getHomogeneousPointsOnSurface(const int N) const {
+  std::vector<Vector> PointsOnSurface;
+  ASSERT(false, "Cuboid_impl::getHomogeneousPointsOnSurface() not implemented yet");
+  return PointsOnSurface;
+}
+
 Shape Cuboid(){
   Shape::impl_ptr impl = Shape::impl_ptr(new Cuboid_impl());
Index: src/Shapes/BaseShapes_impl.hpp
===================================================================
--- src/Shapes/BaseShapes_impl.hpp	(revision 0d02fbf6b0433481fd17dc075f5b57ba0c694fe6)
+++ src/Shapes/BaseShapes_impl.hpp	(revision c5186e909d411ff3c5b6bc38edd607738818eab0)
@@ -13,8 +13,10 @@
 class Sphere_impl : public Shape_impl {
   virtual bool isInside(const Vector &point);
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const;
 };
 
 class Cuboid_impl : public Shape_impl {
   virtual bool isInside(const Vector &point);
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const;
 };
 
Index: src/Shapes/Shape.cpp
===================================================================
--- src/Shapes/Shape.cpp	(revision 0d02fbf6b0433481fd17dc075f5b57ba0c694fe6)
+++ src/Shapes/Shape.cpp	(revision c5186e909d411ff3c5b6bc38edd607738818eab0)
@@ -9,4 +9,8 @@
 #include "Shape_impl.hpp"
 
+
+#include "Helpers/Assert.hpp"
+#include "LinearAlgebra/Vector.hpp"
+
 Shape::Shape(const Shape& src) :
   impl(src.getImpl())
@@ -17,4 +21,8 @@
 bool Shape::isInside(const Vector &point) const{
   return impl->isInside(point);
+}
+
+std::vector<Vector> Shape::getHomogeneousPointsOnSurface(const int N) const {
+  return impl->getHomogeneousPointsOnSurface(N);
 }
 
@@ -65,4 +73,11 @@
 }
 
+std::vector<Vector> AndShape_impl::getHomogeneousPointsOnSurface(const int N) const {
+  std::vector<Vector> PointsOnSurface;
+  ASSERT(false, "AndShape_impl::getHomogeneousPointsOnSurface() not implemented yet");
+  return PointsOnSurface;
+}
+
+
 Shape operator&&(const Shape &lhs,const Shape &rhs){
   Shape::impl_ptr newImpl = Shape::impl_ptr(new AndShape_impl(getShapeImpl(lhs),getShapeImpl(rhs)));
@@ -80,4 +95,21 @@
 bool OrShape_impl::isInside(const Vector &point){
   return rhs->isInside(point) || lhs->isInside(point);
+}
+
+std::vector<Vector> OrShape_impl::getHomogeneousPointsOnSurface(const int N) const {
+  std::vector<Vector> PointsOnSurface_lhs = lhs->getHomogeneousPointsOnSurface(N);
+  std::vector<Vector> PointsOnSurface_rhs = rhs->getHomogeneousPointsOnSurface(N);
+  std::vector<Vector> PointsOnSurface;
+
+  for (std::vector<Vector>::const_iterator iter = PointsOnSurface_lhs.begin(); iter != PointsOnSurface_lhs.end(); ++iter) {
+    if (!rhs->isInside(*iter))
+      PointsOnSurface.push_back(*iter);
+  }
+  for (std::vector<Vector>::const_iterator iter = PointsOnSurface_rhs.begin(); iter != PointsOnSurface_rhs.end(); ++iter) {
+    if (!lhs->isInside(*iter))
+      PointsOnSurface.push_back(*iter);
+  }
+
+  return PointsOnSurface;
 }
 
@@ -99,4 +131,9 @@
 }
 
+std::vector<Vector> NotShape_impl::getHomogeneousPointsOnSurface(const int N) const {
+  // surfaces are the same, only normal direction is different
+  return arg->getHomogeneousPointsOnSurface(N);
+}
+
 Shape operator!(const Shape &arg){
   Shape::impl_ptr newImpl = Shape::impl_ptr(new NotShape_impl(getShapeImpl(arg)));
Index: src/Shapes/Shape.hpp
===================================================================
--- src/Shapes/Shape.hpp	(revision 0d02fbf6b0433481fd17dc075f5b57ba0c694fe6)
+++ src/Shapes/Shape.hpp	(revision c5186e909d411ff3c5b6bc38edd607738818eab0)
@@ -10,4 +10,6 @@
 
 #include <boost/shared_ptr.hpp>
+
+#include <vector>
 
 class Vector;
@@ -25,4 +27,5 @@
 
   bool isInside(const Vector &point) const;
+  std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const;
 
   Shape &operator=(const Shape& rhs);
Index: src/Shapes/ShapeOps.cpp
===================================================================
--- src/Shapes/ShapeOps.cpp	(revision 0d02fbf6b0433481fd17dc075f5b57ba0c694fe6)
+++ src/Shapes/ShapeOps.cpp	(revision c5186e909d411ff3c5b6bc38edd607738818eab0)
@@ -25,4 +25,13 @@
 }
 
+std::vector<Vector> Resize_impl::getHomogeneousPointsOnSurface(const int N) const {
+  std::vector<Vector> PointsOnSurface = arg->getHomogeneousPointsOnSurface(N);
+  for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
+    *iter *= size;
+  }
+  return PointsOnSurface;
+}
+
+
 Shape resize(const Shape &arg,double size){
   Shape::impl_ptr impl = Shape::impl_ptr(new Resize_impl(getShapeImpl(arg),size));
@@ -40,4 +49,12 @@
 bool Translate_impl::isInside(const Vector& point){
   return arg->isInside(point-offset);
+}
+
+std::vector<Vector> Translate_impl::getHomogeneousPointsOnSurface(const int N) const {
+  std::vector<Vector> PointsOnSurface = arg->getHomogeneousPointsOnSurface(N);
+  for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
+    *iter += offset;
+  }
+  return PointsOnSurface;
 }
 
@@ -68,4 +85,12 @@
 }
 
+std::vector<Vector> Stretch_impl::getHomogeneousPointsOnSurface(const int N) const {
+  std::vector<Vector> PointsOnSurface = arg->getHomogeneousPointsOnSurface(N);
+  for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
+    (*iter).ScaleAll(reciFactors);
+  }
+  return PointsOnSurface;
+}
+
 Shape stretch(const Shape &arg, const Vector &factors){
   Shape::impl_ptr impl = Shape::impl_ptr(new Stretch_impl(getShapeImpl(arg),factors));
@@ -87,4 +112,12 @@
 }
 
+std::vector<Vector> Transform_impl::getHomogeneousPointsOnSurface(const int N) const {
+  std::vector<Vector> PointsOnSurface = arg->getHomogeneousPointsOnSurface(N);
+  for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
+    *iter = transformation * (*iter);
+  }
+  return PointsOnSurface;
+}
+
 Shape transform(const Shape &arg, const Matrix &transformation){
   Shape::impl_ptr impl = Shape::impl_ptr(new Transform_impl(getShapeImpl(arg),transformation));
Index: src/Shapes/ShapeOps_impl.hpp
===================================================================
--- src/Shapes/ShapeOps_impl.hpp	(revision 0d02fbf6b0433481fd17dc075f5b57ba0c694fe6)
+++ src/Shapes/ShapeOps_impl.hpp	(revision c5186e909d411ff3c5b6bc38edd607738818eab0)
@@ -19,4 +19,5 @@
   virtual ~Resize_impl();
   virtual bool isInside(const Vector& point);
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const;
 private:
   Shape::impl_ptr arg;
@@ -30,4 +31,5 @@
   virtual ~Translate_impl();
   virtual bool isInside(const Vector& point);
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const;
 private:
   Shape::impl_ptr arg;
@@ -41,4 +43,5 @@
   virtual ~Stretch_impl();
   virtual bool isInside(const Vector& point);
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const;
 private:
   Shape::impl_ptr arg;
@@ -53,4 +56,5 @@
   virtual ~Transform_impl();
   virtual bool isInside(const Vector& point);
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const;
 private:
   Shape::impl_ptr arg;
Index: src/Shapes/Shape_impl.hpp
===================================================================
--- src/Shapes/Shape_impl.hpp	(revision 0d02fbf6b0433481fd17dc075f5b57ba0c694fe6)
+++ src/Shapes/Shape_impl.hpp	(revision c5186e909d411ff3c5b6bc38edd607738818eab0)
@@ -9,5 +9,9 @@
 #define SHAPE_IMPL_HPP_
 
+#include <vector>
+
 #include "Shapes/Shape.hpp"
+
+class Vector;
 
 class Shape_impl {
@@ -16,4 +20,5 @@
   virtual ~Shape_impl(){};
   virtual bool isInside(const Vector &point)=0;
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const=0;
 };
 
@@ -23,4 +28,8 @@
     return true;
   }
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const {
+    std::vector<Vector> PointsOnSurface;
+    return PointsOnSurface;
+  }
 };
 
@@ -28,4 +37,8 @@
   virtual bool isInside(const Vector &point){
     return false;
+  }
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const {
+    std::vector<Vector> PointsOnSurface;
+    return PointsOnSurface;
   }
 };
@@ -36,4 +49,5 @@
   virtual ~AndShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const;
 private:
   Shape::impl_ptr lhs;
@@ -46,4 +60,5 @@
   virtual ~OrShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const;
 private:
   Shape::impl_ptr lhs;
@@ -56,4 +71,5 @@
   virtual ~NotShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const int N) const;
 private:
   Shape::impl_ptr arg;
