Index: src/Shapes/BaseShapes.cpp
===================================================================
--- src/Shapes/BaseShapes.cpp	(revision 9c1c89d06cf8e6779213c5c5181524247fd05e85)
+++ src/Shapes/BaseShapes.cpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -22,9 +22,17 @@
 #include "Shapes/BaseShapes.hpp"
 #include "Shapes/BaseShapes_impl.hpp"
-
-#include <cmath>
+#include "Shapes/ShapeOps.hpp"
+
+#include "defs.hpp"
 
 #include "Helpers/Assert.hpp"
 #include "LinearAlgebra/Vector.hpp"
+#include "LinearAlgebra/Line.hpp"
+#include "LinearAlgebra/Plane.hpp"
+#include "LinearAlgebra/LineSegment.hpp"
+#include "LinearAlgebra/LineSegmentSet.hpp"
+
+#include <cmath>
+#include <algorithm>
 
 bool Sphere_impl::isInside(const Vector &point){
@@ -32,4 +40,27 @@
 }
 
+bool Sphere_impl::isOnSurface(const Vector &point){
+  return fabs(point.NormSquared()-1)<MYEPSILON;
+}
+
+Vector Sphere_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  return point;
+}
+
+LineSegmentSet Sphere_impl::getLineIntersections(const Line &line){
+  LineSegmentSet res(line);
+  std::vector<Vector> intersections = line.getSphereIntersections();
+  if(intersections.size()==2){
+    res.insert(LineSegment(intersections[0],intersections[1]));
+  }
+  return res;
+}
+
+string Sphere_impl::toString(){
+  return "Sphere()";
+}
 
 /**
@@ -66,6 +97,77 @@
 }
 
+Shape Sphere(const Vector &center,double radius){
+  return translate(resize(Sphere(),radius),center);
+}
+
+Shape Ellipsoid(const Vector &center, const Vector &radius){
+  return translate(stretch(Sphere(),radius),center);
+}
+
 bool Cuboid_impl::isInside(const Vector &point){
   return (point[0]>=0 && point[0]<=1) && (point[1]>=0 && point[1]<=1) && (point[2]>=0 && point[2]<=1);
+}
+
+bool Cuboid_impl::isOnSurface(const Vector &point){
+  bool retVal = isInside(point);
+  // test all borders of the cuboid
+  // double fabs
+  retVal = retVal &&
+           (((fabs(point[0]-1.)  < MYEPSILON) || (fabs(point[0])  < MYEPSILON)) ||
+            ((fabs(point[1]-1.)  < MYEPSILON) || (fabs(point[1])  < MYEPSILON)) ||
+            ((fabs(point[2]-1.)  < MYEPSILON) || (fabs(point[2])  < MYEPSILON)));
+  return retVal;
+}
+
+Vector Cuboid_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  Vector res;
+  // figure out on which sides the Vector lies (maximum 3, when it is in a corner)
+  for(int i=NDIM;i--;){
+    if(fabs(fabs(point[i])-1)<MYEPSILON){
+      // add the scaled (-1/+1) Vector to the set of surface vectors
+      res[i] = point[i];
+    }
+  }
+  ASSERT(res.NormSquared()>=1 && res.NormSquared()<=3,"To many or to few sides found for this Vector");
+
+  res.Normalize();
+  return res;
+}
+
+LineSegmentSet Cuboid_impl::getLineIntersections(const Line &line){
+  LineSegmentSet res(line);
+  // get the intersection on each of the six faces
+  vector<Vector> intersections;
+  intersections.resize(2);
+  int c=0;
+  int x[2]={-1,+1};
+  for(int i=NDIM;i--;){
+    for(int p=0;p<2;++p){
+      if(c==2) goto end; // I know this sucks, but breaking two loops is stupid
+      Vector base;
+      base[i]=x[p];
+      // base now points to the surface and is normal to it at the same time
+      Plane p(base,base);
+      Vector intersection = p.GetIntersection(line);
+      if(isInside(intersection)){
+        // if we have a point on the edge it might already be contained
+        if(c==1 && intersections[0]==intersection)
+          continue;
+        intersections[c++]=intersection;
+      }
+    }
+  }
+  end:
+  if(c==2){
+    res.insert(LineSegment(intersections[0],intersections[1]));
+  }
+  return res;
+}
+
+string Cuboid_impl::toString(){
+  return "Cuboid()";
 }
 
@@ -83,2 +185,17 @@
   return Shape(impl);
 }
+
+Shape Cuboid(const Vector &corner1, const Vector &corner2){
+  // make sure the two edges are upper left front and lower right back
+  Vector sortedC1;
+  Vector sortedC2;
+  for(int i=NDIM;i--;){
+    sortedC1[i] = min(corner1[i],corner2[i]);
+    sortedC2[i] = max(corner1[i],corner2[i]);
+    ASSERT(corner1[i]!=corner2[i],"Given points for cuboid edges did not define a valid space");
+  }
+  // get the middle point
+  Vector middle = (1./2.)*(sortedC1+sortedC2);
+  Vector factors = sortedC2-middle;
+  return translate(stretch(Cuboid(),factors),middle);
+}
Index: src/Shapes/BaseShapes.hpp
===================================================================
--- src/Shapes/BaseShapes.hpp	(revision 9c1c89d06cf8e6779213c5c5181524247fd05e85)
+++ src/Shapes/BaseShapes.hpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -12,5 +12,8 @@
 
 Shape Sphere();
+Shape Sphere(const Vector &center,double radius);
+Shape Ellipsoid(const Vector &center, const Vector &radius);
 Shape Cuboid();
+Shape Cuboid(const Vector &corner1, const Vector &corner2);
 
 #endif /* BASESHAPES_HPP_ */
Index: src/Shapes/BaseShapes_impl.hpp
===================================================================
--- src/Shapes/BaseShapes_impl.hpp	(revision 9c1c89d06cf8e6779213c5c5181524247fd05e85)
+++ src/Shapes/BaseShapes_impl.hpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -9,8 +9,15 @@
 #define BASESHAPES_IMPL_HPP_
 
+class Vector;
+class LineSegmentSet;
+
 #include "Shapes/Shape_impl.hpp"
 
 class Sphere_impl : public Shape_impl {
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 };
@@ -18,4 +25,8 @@
 class Cuboid_impl : public Shape_impl {
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 };
Index: src/Shapes/Shape.cpp
===================================================================
--- src/Shapes/Shape.cpp	(revision 9c1c89d06cf8e6779213c5c5181524247fd05e85)
+++ src/Shapes/Shape.cpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -20,10 +20,12 @@
 #include "Helpers/MemDebug.hpp"
 
+#include "Helpers/Assert.hpp"
+#include "LinearAlgebra/Vector.hpp"
+
 #include "Shape.hpp"
 #include "Shape_impl.hpp"
 
-
-#include "Helpers/Assert.hpp"
-#include "LinearAlgebra/Vector.hpp"
+#include <string>
+
 
 Shape::Shape(const Shape& src) :
@@ -35,4 +37,16 @@
 bool Shape::isInside(const Vector &point) const{
   return impl->isInside(point);
+}
+
+bool Shape::isOnSurface(const Vector &point) const{
+  return impl->isOnSurface(point);
+}
+
+Vector Shape::getNormal(const Vector &point) const throw (NotOnSurfaceException){
+  return impl->getNormal(point);
+}
+
+LineSegmentSet Shape::getLineIntersections(const Line &line){
+  return impl->getLineIntersections(line);
 }
 
@@ -52,4 +66,8 @@
 }
 
+std::string Shape::toString() const{
+  return impl->toString();
+}
+
 Shape::impl_ptr Shape::getImpl() const{
   return impl;
@@ -85,4 +103,55 @@
 bool AndShape_impl::isInside(const Vector &point){
   return lhs->isInside(point) && rhs->isInside(point);
+}
+
+bool AndShape_impl::isOnSurface(const Vector &point){
+  // check the number of surfaces that this point is on
+  int surfaces =0;
+  surfaces += lhs->isOnSurface(point);
+  surfaces += rhs->isOnSurface(point);
+
+  switch(surfaces){
+    case 0:
+      return false;
+      // no break necessary
+    case 1:
+      // if it is inside for the object where it does not lie on
+      // the surface the whole point lies inside
+      return (lhs->isOnSurface(point) && rhs->isInside(point)) ||
+             (rhs->isOnSurface(point) && lhs->isInside(point));
+      // no break necessary
+    case 2:
+      {
+        // it lies on both Shapes... could be an edge or an inner point
+        // test the direction of the normals
+        Vector direction=lhs->getNormal(point)+rhs->getNormal(point);
+        // if the directions are opposite we lie on the inside
+        return !direction.IsZero();
+      }
+      // no break necessary
+    default:
+      // if this happens there is something wrong
+      ASSERT(0,"Default case should have never been used");
+  }
+  return false; // never reached
+}
+
+Vector AndShape_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
+  Vector res;
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  res += lhs->isOnSurface(point)?lhs->getNormal(point):zeroVec;
+  res += rhs->isOnSurface(point)?rhs->getNormal(point):zeroVec;
+  res.Normalize();
+  return res;
+}
+
+LineSegmentSet AndShape_impl::getLineIntersections(const Line &line){
+  return intersect(lhs->getLineIntersections(line),rhs->getLineIntersections(line));
+}
+
+string AndShape_impl::toString(){
+  return string("(") + lhs->toString() + string("&&") + rhs->toString() + string(")");
 }
 
@@ -122,4 +191,55 @@
 }
 
+bool OrShape_impl::isOnSurface(const Vector &point){
+  // check the number of surfaces that this point is on
+  int surfaces =0;
+  surfaces += lhs->isOnSurface(point);
+  surfaces += rhs->isOnSurface(point);
+
+  switch(surfaces){
+    case 0:
+      return false;
+      // no break necessary
+    case 1:
+      // if it is inside for the object where it does not lie on
+      // the surface the whole point lies inside
+      return (lhs->isOnSurface(point) && !rhs->isInside(point)) ||
+             (rhs->isOnSurface(point) && !lhs->isInside(point));
+      // no break necessary
+    case 2:
+      {
+        // it lies on both Shapes... could be an edge or an inner point
+        // test the direction of the normals
+        Vector direction=lhs->getNormal(point)+rhs->getNormal(point);
+        // if the directions are opposite we lie on the inside
+        return !direction.IsZero();
+      }
+      // no break necessary
+    default:
+      // if this happens there is something wrong
+      ASSERT(0,"Default case should have never been used");
+  }
+  return false; // never reached
+}
+
+Vector OrShape_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
+  Vector res;
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  res += lhs->isOnSurface(point)?lhs->getNormal(point):zeroVec;
+  res += rhs->isOnSurface(point)?rhs->getNormal(point):zeroVec;
+  res.Normalize();
+  return res;
+}
+
+LineSegmentSet OrShape_impl::getLineIntersections(const Line &line){
+  return merge(lhs->getLineIntersections(line),rhs->getLineIntersections(line));
+}
+
+string OrShape_impl::toString(){
+  return string("(") + lhs->toString() + string("||") + rhs->toString() + string(")");
+}
+
 std::vector<Vector> OrShape_impl::getHomogeneousPointsOnSurface(const size_t N) const {
   std::vector<Vector> PointsOnSurface_lhs = lhs->getHomogeneousPointsOnSurface(N);
@@ -156,4 +276,20 @@
 }
 
+bool NotShape_impl::isOnSurface(const Vector &point){
+  return arg->isOnSurface(point);
+}
+
+Vector NotShape_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
+  return -1*arg->getNormal(point);
+}
+
+LineSegmentSet NotShape_impl::getLineIntersections(const Line &line){
+  return invert(arg->getLineIntersections(line));
+}
+
+string NotShape_impl::toString(){
+  return string("!") + arg->toString();
+}
+
 std::vector<Vector> NotShape_impl::getHomogeneousPointsOnSurface(const size_t N) const {
   // surfaces are the same, only normal direction is different
@@ -165,2 +301,8 @@
   return Shape(newImpl);
 }
+
+/**************** global operations *********************************/
+ostream &operator<<(ostream &ost,const Shape &shape){
+  ost << shape.toString();
+  return ost;
+}
Index: src/Shapes/Shape.hpp
===================================================================
--- src/Shapes/Shape.hpp	(revision 9c1c89d06cf8e6779213c5c5181524247fd05e85)
+++ src/Shapes/Shape.hpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -10,4 +10,7 @@
 
 #include <boost/shared_ptr.hpp>
+#include <iosfwd>
+
+#include "Exceptions/NotOnSurfaceException.hpp"
 
 #include <vector>
@@ -15,4 +18,6 @@
 class Vector;
 class Shape_impl;
+class LineSegmentSet;
+class Line;
 
 class Shape
@@ -27,8 +32,13 @@
 
   bool isInside(const Vector &point) const;
+  bool isOnSurface(const Vector &point) const;
+  Vector getNormal(const Vector &point) const throw(NotOnSurfaceException);
+
+  LineSegmentSet getLineIntersections(const Line&);
   std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 
   Shape &operator=(const Shape& rhs);
 
+  std::string toString() const;
 protected:
   impl_ptr getImpl() const;
@@ -45,3 +55,5 @@
 Shape operator!(const Shape&);
 
+std::ostream &operator<<(std::ostream&,const Shape&);
+
 #endif /* SHAPE_HPP_ */
Index: src/Shapes/ShapeOps.cpp
===================================================================
--- src/Shapes/ShapeOps.cpp	(revision 9c1c89d06cf8e6779213c5c5181524247fd05e85)
+++ src/Shapes/ShapeOps.cpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -23,10 +23,63 @@
 #include "Shapes/ShapeOps_impl.hpp"
 
+#include "LinearAlgebra/Vector.hpp"
 #include "Helpers/Assert.hpp"
 
+/*************** Base case ***********************/
+
+ShapeOpsBase_impl::ShapeOpsBase_impl(const Shape::impl_ptr &_arg) :
+  arg(_arg){}
+
+ShapeOpsBase_impl::~ShapeOpsBase_impl(){}
+
+bool ShapeOpsBase_impl::isInside(const Vector &point){
+  return arg->isInside(translateIn(point));
+}
+
+bool ShapeOpsBase_impl::isOnSurface(const Vector &point){
+  return arg->isOnSurface(translateIn(point));
+}
+
+Vector ShapeOpsBase_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
+  Vector helper = translateIn(point);
+  if(!arg->isOnSurface(helper)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  Vector res = translateOutNormal(arg->getNormal(helper));
+  res.Normalize();
+  return res;
+}
+
+LineSegmentSet ShapeOpsBase_impl::getLineIntersections(const Line &line){
+  Line newLine(translateIn(line.getOrigin()),translateIn(line.getDirection()));
+  LineSegmentSet res(line);
+  LineSegmentSet helper = getArg()->getLineIntersections(newLine);
+  for(LineSegmentSet::iterator iter = helper.begin();iter!=helper.end();++iter){
+    LinePoint lpBegin = iter->getBegin();
+    LinePoint lpEnd = iter->getBegin();
+    // translate both linepoints
+    lpBegin = lpBegin.isNegInfinity()?
+                line.negEndpoint():
+                line.getLinePoint(translateOutPos(lpBegin.getPoint()));
+    lpEnd = lpEnd.isPosInfinity()?
+              line.posEndpoint():
+              line.getLinePoint(translateOutPos(lpEnd.getPoint()));
+    res.insert(LineSegment(lpBegin,lpEnd));
+  }
+  return res;
+}
+
+std::vector<Vector> ShapeOpsBase_impl::getHomogeneousPointsOnSurface(const size_t N) const {
+  return getArg()->getHomogeneousPointsOnSurface(N);;
+}
+
+Shape::impl_ptr ShapeOpsBase_impl::getArg() const{
+  return arg;
+}
+
 /********************* Resize ********************/
 
 Resize_impl::Resize_impl(const Shape::impl_ptr &_arg,double _size) :
-  arg(_arg), size(_size)
+  ShapeOpsBase_impl(_arg), size(_size)
 {
   ASSERT(size>0,"Cannot resize a Shape to size zero or below");
@@ -36,9 +89,27 @@
 
 bool Resize_impl::isInside(const Vector& point){
-  return arg->isInside((1/size) * point);
+  return getArg()->isInside((1/size) * point);
+}
+
+Vector Resize_impl::translateIn(const Vector& point){
+  return (1/size) * point;
+}
+
+Vector Resize_impl::translateOutPos(const Vector& point){
+  return size * point;
+}
+
+Vector Resize_impl::translateOutNormal(const Vector& point){
+  return point;
+}
+
+string Resize_impl::toString(){
+  stringstream sstr;
+  sstr << "resize(" << getArg()->toString() << "," << size << ")";
+  return sstr.str();
 }
 
 std::vector<Vector> Resize_impl::getHomogeneousPointsOnSurface(const size_t N) const {
-  std::vector<Vector> PointsOnSurface = arg->getHomogeneousPointsOnSurface(N);
+  std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
   for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
     *iter *= size;
@@ -56,5 +127,5 @@
 
 Translate_impl::Translate_impl(const Shape::impl_ptr &_arg, const Vector &_offset) :
-  arg(_arg),offset(_offset)
+  ShapeOpsBase_impl(_arg),offset(_offset)
 {}
 
@@ -62,9 +133,27 @@
 
 bool Translate_impl::isInside(const Vector& point){
-  return arg->isInside(point-offset);
+  return getArg()->isInside(point-offset);
+}
+
+Vector Translate_impl::translateIn(const Vector& point){
+  return point-offset;
+}
+
+Vector Translate_impl::translateOutPos(const Vector& point){
+  return point+offset;
+}
+
+Vector Translate_impl::translateOutNormal(const Vector& point){
+  return point;
+}
+
+string Translate_impl::toString(){
+  stringstream sstr;
+  sstr << "translate(" << getArg()->toString() << "," << offset << ")";
+  return sstr.str();
 }
 
 std::vector<Vector> Translate_impl::getHomogeneousPointsOnSurface(const size_t N) const {
-  std::vector<Vector> PointsOnSurface = arg->getHomogeneousPointsOnSurface(N);
+  std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
   for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
     *iter += offset;
@@ -81,5 +170,5 @@
 
 Stretch_impl::Stretch_impl(const Shape::impl_ptr &_arg, const Vector &_factors) :
-  arg(_arg),factors(_factors)
+  ShapeOpsBase_impl(_arg),factors(_factors)
 {
   ASSERT(factors[0]>0,"cannot stretch a shape by a negative amount");
@@ -96,9 +185,39 @@
   Vector helper=point;
   helper.ScaleAll(reciFactors);
-  return arg->isInside(helper);
+  return getArg()->isInside(helper);
+}
+
+Vector Stretch_impl::translateIn(const Vector& point){
+  Vector helper=point;
+  helper.ScaleAll(reciFactors);
+  return helper;
+}
+
+Vector Stretch_impl::translateOutPos(const Vector& point){
+  Vector helper=point;
+  helper.ScaleAll(factors);
+  return helper;
+}
+
+Vector Stretch_impl::translateOutNormal(const Vector& point){
+  Vector helper=point;
+  // the normalFactors are derived from appearances of the factors
+  // with in the vectorproduct
+  Vector normalFactors;
+  normalFactors[0]=factors[1]*factors[2];
+  normalFactors[1]=factors[0]*factors[2];
+  normalFactors[2]=factors[0]*factors[1];
+  helper.ScaleAll(normalFactors);
+  return helper;
+}
+
+string Stretch_impl::toString(){
+  stringstream sstr;
+  sstr << "stretch(" << getArg()->toString() << "," << factors << ")";
+  return sstr.str();
 }
 
 std::vector<Vector> Stretch_impl::getHomogeneousPointsOnSurface(const size_t N) const {
-  std::vector<Vector> PointsOnSurface = arg->getHomogeneousPointsOnSurface(N);
+  std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
   for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
     (*iter).ScaleAll(reciFactors);
@@ -115,5 +234,5 @@
 
 Transform_impl::Transform_impl(const Shape::impl_ptr &_arg, const Matrix &_transformation) :
-  arg(_arg),transformation(_transformation)
+  ShapeOpsBase_impl(_arg),transformation(_transformation)
 {
   transformationInv = transformation.invert();
@@ -123,9 +242,28 @@
 
 bool Transform_impl::isInside(const Vector& point){
-  return arg->isInside(transformationInv * point);
+  return getArg()->isInside(transformationInv * point);
+}
+
+Vector Transform_impl::translateIn(const Vector& point){
+  return transformationInv * point;
+}
+
+Vector Transform_impl::translateOutPos(const Vector& point){
+  return transformation * point;
+}
+
+Vector Transform_impl::translateOutNormal(const Vector& point){
+  Matrix mat = transformation.invert().transpose();
+  return mat * point;
+}
+
+string Transform_impl::toString(){
+  stringstream sstr;
+  sstr << "transform(" << getArg()->toString() << "," << transformation << ")";
+  return sstr.str();
 }
 
 std::vector<Vector> Transform_impl::getHomogeneousPointsOnSurface(const size_t N) const {
-  std::vector<Vector> PointsOnSurface = arg->getHomogeneousPointsOnSurface(N);
+  std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
   for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
     *iter = transformation * (*iter);
Index: src/Shapes/ShapeOps_impl.hpp
===================================================================
--- src/Shapes/ShapeOps_impl.hpp	(revision 9c1c89d06cf8e6779213c5c5181524247fd05e85)
+++ src/Shapes/ShapeOps_impl.hpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -13,50 +13,88 @@
 #include "LinearAlgebra/Matrix.hpp"
 
-class Resize_impl :  public Shape_impl
+#include <vector>
+
+class LineSegment;
+
+class ShapeOpsBase_impl : public Shape_impl{
+public:
+  ShapeOpsBase_impl(const Shape::impl_ptr&);
+  virtual ~ShapeOpsBase_impl();
+  virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw (NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
+protected:
+  virtual Vector translateIn(const Vector &point)=0;
+  virtual Vector translateOutPos(const Vector &point)=0;
+  virtual Vector translateOutNormal(const Vector &point)=0;
+  Shape::impl_ptr getArg() const;
+private:
+  Shape::impl_ptr arg;
+};
+
+class Resize_impl :  public ShapeOpsBase_impl
 {
 public:
   Resize_impl(const Shape::impl_ptr&,double);
   virtual ~Resize_impl();
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
+  virtual std::string toString();
   virtual bool isInside(const Vector& point);
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
-  Shape::impl_ptr arg;
   double size;
 };
 
-class Translate_impl :  public Shape_impl
+class Translate_impl :  public ShapeOpsBase_impl
 {
 public:
   Translate_impl(const Shape::impl_ptr&, const Vector&);
   virtual ~Translate_impl();
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
+  virtual std::string toString();
   virtual bool isInside(const Vector& point);
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
-  Shape::impl_ptr arg;
   Vector offset;
 };
 
-class Stretch_impl : public Shape_impl
+class Stretch_impl : public ShapeOpsBase_impl
 {
 public:
   Stretch_impl(const Shape::impl_ptr&, const Vector&);
   virtual ~Stretch_impl();
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
+  virtual std::string toString();
   virtual bool isInside(const Vector& point);
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
-  Shape::impl_ptr arg;
   Vector factors;
   Vector reciFactors;
 };
 
-class Transform_impl : public Shape_impl
+class Transform_impl : public ShapeOpsBase_impl
 {
 public:
   Transform_impl(const Shape::impl_ptr&, const Matrix&);
   virtual ~Transform_impl();
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
+  virtual std::string toString();
   virtual bool isInside(const Vector& point);
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
-  Shape::impl_ptr arg;
   Matrix transformation;
   Matrix transformationInv;
Index: src/Shapes/Shape_impl.hpp
===================================================================
--- src/Shapes/Shape_impl.hpp	(revision 9c1c89d06cf8e6779213c5c5181524247fd05e85)
+++ src/Shapes/Shape_impl.hpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -12,6 +12,9 @@
 
 #include "Shapes/Shape.hpp"
+#include "LinearAlgebra/Line.hpp"
+#include "LinearAlgebra/LineSegment.hpp"
+#include "LinearAlgebra/LineSegmentSet.hpp"
+#include "LinearAlgebra/Vector.hpp"
 
-class Vector;
 
 class Shape_impl {
@@ -20,4 +23,8 @@
   virtual ~Shape_impl(){};
   virtual bool isInside(const Vector &point)=0;
+  virtual bool isOnSurface(const Vector &point)=0;
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException)=0;
+  virtual LineSegmentSet getLineIntersections(const Line&)=0;
+  virtual std::string toString()=0;
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const=0;
 };
@@ -27,4 +34,18 @@
   virtual bool isInside(const Vector &point){
     return true;
+  }
+  virtual bool isOnSurface(const Vector &point){
+    return false;
+  }
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  virtual LineSegmentSet getLineIntersections(const Line &line){
+    LineSegmentSet res(line);
+    res.insert(LineSegment(line.negEndpoint(),line.posEndpoint()));
+    return res;
+  }
+  virtual std::string toString(){
+    return "Everywhere()";
   }
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const {
@@ -37,4 +58,16 @@
   virtual bool isInside(const Vector &point){
     return false;
+  }
+  virtual bool isOnSurface(const Vector &point){
+    return false;
+  }
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  virtual LineSegmentSet getLineIntersections(const Line &line){
+    return LineSegmentSet(line);
+  }
+  virtual std::string toString(){
+    return "Nowhere()";
   }
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const {
@@ -49,4 +82,8 @@
   virtual ~AndShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
@@ -60,4 +97,8 @@
   virtual ~OrShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
@@ -71,4 +112,8 @@
   virtual ~NotShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
   virtual std::vector<Vector> getHomogeneousPointsOnSurface(const size_t N) const;
 private:
