soc-2008-mxcurioni: merged changes to revision 14747, cosmetic changes for source/blender/freestyle

27 files changed, 6086 insertions, 0 deletions

diff --git a/source/blender/freestyle/intern/geometry/BBox.h b/source/blender/freestyle/intern/geometry/BBox.h
new file mode 100755
index 00000000000..9c46d7918e2
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/BBox.h
@@ -0,0 +1,141 @@
+//
+//  Filename         : BBox.h
+//  Author(s)        : Stephane Grabli
+//  Purpose          : A class to hold a bounding box
+//  Date of creation : 22/05/2003
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef  BBOX_H
+# define BBOX_H
+
+template <class Point>
+class BBox
+{
+ public:
+
+  inline BBox() {
+    _empty = true;
+  }
+
+  template <class T>
+  inline BBox(const T& min_in, const T& max_in) : _min(min_in), _max(max_in) {
+    _empty = false;
+  }
+
+  template <class T>
+  inline BBox(const BBox<T>& b) : _min(b.getMin()), _max(b.getMax()) {
+    _empty = false;
+  }
+
+  template <class T>
+  inline void extendToContain(const T& p) {
+    if (_empty) {
+      _min = p;
+      _max = p;
+      _empty = false;
+      return;
+    }
+    for (unsigned i = 0; i < Point::dim(); i++) {
+      if (p[i] < _min[i])
+	_min[i] = p[i];
+      else if (p[i] > _max[i])
+	_max[i] = p[i];
+    }
+    _empty = false;
+  }
+
+  inline void clear() {
+    _empty = true;
+  }
+
+  inline bool empty() const {
+    return _empty;
+  }
+
+  inline const Point& getMin() const {
+    return _min;
+  }
+
+  inline const Point& getMax() const {
+    return _max;
+  }
+
+  inline BBox<Point>& operator=(const BBox<Point>& b) {
+    _min = b.getMin();
+    _max = b.getMax();
+    _empty = false;
+    return *this;
+  }
+
+  inline BBox<Point>& operator+=(const BBox<Point>& b) {
+    if (_empty) {
+      _min = b.getMin();
+      _max = b.getMax();
+      _empty = false;
+    }
+    else {
+      for (unsigned i = 0; i < Point::dim(); i++) {
+	if (b.getMin()[i] < _min[i])
+	  _min[i] = b.getMin()[i];
+	if (b.getMax()[i] > _max[i])
+	  _max[i] = b.getMax()[i];
+      }
+    }
+    return *this;
+  }
+
+  inline bool inside(const Point& p){
+    if(empty())
+        return false;
+    for (unsigned i = 0; i < Point::dim(); i++) {
+        if((_min[i]>p[i]) || (_max[i]<p[i]))
+            return false;
+    }
+    return true;
+
+  }
+
+private:
+
+  Point	_min;
+  Point	_max;
+  bool	_empty;
+};
+
+template <class Point>
+BBox<Point>& operator+(const BBox<Point> &b1, const BBox<Point> &b2)
+{
+  Point new_min;
+  Point new_max;
+
+  for (unsigned i = 0; i < Point::dim(); i++) {
+    new_min[i] = b1.getMin()[i] < b2.getMin()[i] ? b1.getMin()[i] : b2.getMin()[i];
+    new_max[i] = b1.getMax()[i] > b2.getMax()[i] ? b1.getMax()[i] : b2.getMax()[i];
+  }
+
+  return BBox<Point>(new_min, new_max);
+}
+
+#endif // BBOX_H
diff --git a/source/blender/freestyle/intern/geometry/Bezier.cpp b/source/blender/freestyle/intern/geometry/Bezier.cpp
new file mode 100755
index 00000000000..8f9771f29d3
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/Bezier.cpp
@@ -0,0 +1,118 @@
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#include "Bezier.h"
+#include "FitCurve.h"
+
+using namespace std;
+
+BezierCurveSegment::BezierCurveSegment()
+{
+}
+
+BezierCurveSegment::~BezierCurveSegment()
+{
+}
+
+void BezierCurveSegment::AddControlPoint(const Vec2d& iPoint)
+{
+  _ControlPolygon.push_back(iPoint);
+  if(_ControlPolygon.size() == 4)
+    Build();
+}
+
+void BezierCurveSegment::Build()
+{
+  if(_ControlPolygon.size() != 4)
+    return;
+
+  // Compute the rightmost part of the matrix:
+  vector<Vec2d>::const_iterator p0,p1,p2,p3;
+  p0 = _ControlPolygon.begin();
+  p1 = p0;++p1;
+  p2 = p1;++p2;
+  p3 = p2;++p3;
+  float x[4], y[4];
+  
+  x[0] = -p0->x()+3*p1->x()-3*p2->x()+p3->x();
+  x[1] = 3*p0->x()-6*p1->x()+3*p2->x();
+  x[2] = -3*p0->x()+3*p1->x();
+  x[3] = p0->x();
+
+  y[0] = -p0->y()+3*p1->y()-3*p2->y()+p3->y();
+  y[1] = 3*p0->y()-6*p1->y()+3*p2->y();
+  y[2] = -3*p0->y()+3*p1->y();
+  y[3] = p0->y();
+
+  int nvertices = 12;
+  float increment = 1.0/(float)nvertices;
+  float t = 0.f;
+  for(int i=0; i<=nvertices; ++i)
+  {
+    _Vertices.push_back(Vec2d((x[3] + t*(x[2] + t*(x[1] + t*x[0]))), 
+                              (y[3] + t*(y[2] + t*(y[1] + t*y[0])))));
+    t+=increment;
+  }
+}
+
+BezierCurve::BezierCurve()
+{
+  _currentSegment = new BezierCurveSegment;
+}
+
+BezierCurve::BezierCurve(vector<Vec2d>& iPoints, double error)
+{
+  FitCurveWrapper fitcurve;
+  _currentSegment = new BezierCurveSegment;
+  vector<Vec2d> curve;
+
+  fitcurve.FitCurve(iPoints, curve, error);
+  int i=0;
+  vector<Vec2d>::iterator v,vend;
+  for(v=curve.begin(),vend=curve.end();
+  v!=vend;
+  ++v)
+  {
+    if((i == 0) || (i%4 != 0))
+      AddControlPoint(*v);
+    ++i;
+  }
+}
+
+BezierCurve::~BezierCurve()
+{
+  if(_currentSegment)
+    delete _currentSegment;
+}
+
+void BezierCurve::AddControlPoint(const Vec2d& iPoint)
+{
+  _ControlPolygon.push_back(iPoint);
+  _currentSegment->AddControlPoint(iPoint);
+  if(_currentSegment->size() == 4)
+  {
+    _Segments.push_back(_currentSegment);
+    _currentSegment = new BezierCurveSegment;
+    _currentSegment->AddControlPoint(iPoint);
+  }
+}
+
+
diff --git a/source/blender/freestyle/intern/geometry/Bezier.h b/source/blender/freestyle/intern/geometry/Bezier.h
new file mode 100755
index 00000000000..acae71bbb2c
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/Bezier.h
@@ -0,0 +1,73 @@
+//
+//  Filename         : Bezier.h
+//  Author(s)        : Stephane Grabli
+//  Purpose          : Class to define a Bezier curve of order 4.
+//  Date of creation : 04/06/2003
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef  BEZIER_H
+# define BEZIER_H
+
+#include <vector>
+#include "../system/FreestyleConfig.h"
+#include "Geom.h"
+
+using namespace Geometry;
+
+class LIB_GEOMETRY_EXPORT BezierCurveSegment 
+{
+private:
+  std::vector<Vec2d> _ControlPolygon;
+  std::vector<Vec2d> _Vertices;
+
+public:
+  BezierCurveSegment();
+  virtual ~BezierCurveSegment();
+
+  void AddControlPoint(const Vec2d& iPoint);
+  void Build();
+  inline int size() const {return _ControlPolygon.size();}
+  inline std::vector<Vec2d>& vertices() {return _Vertices;}
+};
+
+
+class LIB_GEOMETRY_EXPORT BezierCurve 
+{
+private:
+  std::vector<Vec2d> _ControlPolygon;
+  std::vector<BezierCurveSegment*> _Segments;
+  BezierCurveSegment *_currentSegment;
+
+public:
+  BezierCurve();
+  BezierCurve(std::vector<Vec2d>& iPoints, double error=4.0);
+  virtual ~BezierCurve();
+
+  void AddControlPoint(const Vec2d& iPoint);
+  std::vector<Vec2d>& controlPolygon() {return _ControlPolygon;}
+  std::vector<BezierCurveSegment*>& segments() {return _Segments;}
+};
+
+#endif // BEZIER_H
diff --git a/source/blender/freestyle/intern/geometry/FastGrid.cpp b/source/blender/freestyle/intern/geometry/FastGrid.cpp
new file mode 100755
index 00000000000..b090a3df67f
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/FastGrid.cpp
@@ -0,0 +1,62 @@
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#include "FastGrid.h"
+
+void FastGrid::clear() {
+  if(!_cells)
+    return;
+
+  for(unsigned i = 0; i < _cells_size; i++)
+    if (_cells[i])
+      delete _cells[i];
+  delete[] _cells;
+  _cells = NULL;
+  _cells_size = 0;
+
+  Grid::clear();
+}
+
+void FastGrid::configure(const Vec3r& orig, const Vec3r& size, unsigned nb) {
+  Grid::configure(orig, size, nb);
+  _cells_size = _cells_nb[0] * _cells_nb[1] * _cells_nb[2];
+  _cells = new Cell*[_cells_size];
+  memset(_cells, 0, _cells_size * sizeof(*_cells));
+}
+
+Cell* FastGrid::getCell(const Vec3u& p) {
+    //cout << _cells<< " "<< p << " " <<_cells_nb[0]<<"-"<< _cells_nb[1]<<"-"<< _cells_nb[2]<< " "<<_cells_size<< endl;
+    assert(_cells||("_cells is a null pointer"));
+    assert((_cells_nb[0] * (p[2] * _cells_nb[1] + p[1]) + p[0])<_cells_size);
+    assert(p[0]<_cells_nb[0]);
+    assert(p[1]<_cells_nb[1]);
+    assert(p[2]<_cells_nb[2]);
+    return _cells[_cells_nb[0] * (p[2] * _cells_nb[1] + p[1]) + p[0]];
+}
+
+void FastGrid::fillCell(const Vec3u& p, Cell& cell) {
+    assert(_cells||("_cells is a null pointer"));
+    assert((_cells_nb[0] * (p[2] * _cells_nb[1] + p[1]) + p[0])<_cells_size);
+    assert(p[0]<_cells_nb[0]);
+    assert(p[1]<_cells_nb[1]);
+    assert(p[2]<_cells_nb[2]);
+    _cells[_cells_nb[0] * (p[2] * _cells_nb[1] + p[1]) + p[0]] = &cell;
+}
\ No newline at end of file
diff --git a/source/blender/freestyle/intern/geometry/FastGrid.h b/source/blender/freestyle/intern/geometry/FastGrid.h
new file mode 100755
index 00000000000..e620ff24385
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/FastGrid.h
@@ -0,0 +1,85 @@
+//
+//  Filename         : FastGrid.h
+//  Author(s)        : Stephane Grabli
+//  Purpose          : Class to define a cell grid surrounding the
+//                     bounding box of the scene
+//  Date of creation : 30/07/2002
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef  FASTGRID_H
+# define FASTGRID_H
+
+# include "Grid.h"
+# include <cassert>
+/*! Class to define a regular grid used for ray
+ *  casting computations 
+ *  We don't use a hashtable here. The grid is 
+ *  explicitly stored for faster computations.
+ *  However, this might result in significant
+ *  increase in memory usage (compared to the regular grid)
+ */
+
+class LIB_GEOMETRY_EXPORT FastGrid : public Grid
+{
+ public:
+
+  FastGrid() : Grid() {
+    _cells = NULL;
+    _cells_size = 0;
+  }
+  
+  virtual ~FastGrid() {
+    clear();
+  }
+
+  /*! clears the grid
+   *  Deletes all the cells, clears the hashtable,
+   *  resets size, size of cell, number of cells.
+   */
+  virtual void clear();
+  
+  /*! Sets the different parameters of the grid
+   *    orig
+   *      The grid origin
+   *    size
+   *      The grid's dimensions
+   *    nb
+   *      The number of cells of the grid
+   */
+  virtual void configure(const Vec3r& orig, const Vec3r& size, unsigned nb);
+  
+  /*! returns the cell whose coordinates are pased as argument */
+  Cell* getCell(const Vec3u& p) ;
+ 
+  /*! Fills the case p with the cell iCell */
+  virtual void fillCell(const Vec3u& p, Cell& cell);
+
+protected:
+
+  Cell**	_cells;
+  unsigned	_cells_size;
+};
+
+#endif // FASTGRID_H
diff --git a/source/blender/freestyle/intern/geometry/FitCurve.cpp b/source/blender/freestyle/intern/geometry/FitCurve.cpp
new file mode 100755
index 00000000000..ade40b050ca
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/FitCurve.cpp
@@ -0,0 +1,602 @@
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#include <stdio.h>
+#include <math.h>
+#include "FitCurve.h"
+
+using namespace std;
+
+typedef Vector2 *BezierCurve;
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
+/* Forward declarations */
+static double *Reparameterize(Vector2 *d, int first, int last, double *u, BezierCurve bezCurve);
+static double NewtonRaphsonRootFind(BezierCurve Q, Vector2 P, double u);
+static Vector2 BezierII(int degree, Vector2 *V, double t);
+static double B0(double u);
+static double B1(double u);
+static double B2(double u);
+static double B3(double u);
+static Vector2 ComputeLeftTangent(Vector2 *d, int end);
+static Vector2 ComputeLeftTangent(Vector2 *d, int end);
+static Vector2 ComputeLeftTangent(Vector2 *d, int end);
+static double ComputeMaxError(Vector2 *d, int first, int last, BezierCurve bezCurve, double *u, int *splitPoint);
+static double *ChordLengthParameterize(Vector2 *d, int first, int last);
+static BezierCurve  GenerateBezier(Vector2 *d, int first, int last, double *uPrime, Vector2 tHat1, Vector2 tHat2);
+static Vector2 V2AddII(Vector2 a, Vector2 b);
+static Vector2 V2ScaleIII(Vector2 v, double s);
+static Vector2 V2SubII(Vector2 a, Vector2 b);
+
+
+#define MAXPOINTS	1000		/* The most points you can have */
+
+/* returns squared length of input vector */	
+double V2SquaredLength(Vector2 *a) 
+{	return(((*a)[0] * (*a)[0])+((*a)[1] * (*a)[1]));
+}
+	
+/* returns length of input vector */
+double V2Length(Vector2 *a) 
+{
+	return(sqrt(V2SquaredLength(a)));
+}
+
+Vector2 *V2Scale(Vector2 *v, double newlen) 
+{
+  double len = V2Length(v);
+	if (len != 0.0) { (*v)[0] *= newlen/len;   (*v)[1] *= newlen/len; }
+	return(v);
+}
+
+/* return the dot product of vectors a and b */
+double V2Dot(Vector2 *a, Vector2 *b) 
+{
+	return(((*a)[0]*(*b)[0])+((*a)[1]*(*b)[1]));
+}
+
+/* return the distance between two points */
+double V2DistanceBetween2Points(Vector2 *a, Vector2 *b)
+{
+double dx = (*a)[0] - (*b)[0];
+double dy = (*a)[1] - (*b)[1];
+	return(sqrt((dx*dx)+(dy*dy)));
+}
+
+/* return vector sum c = a+b */
+Vector2 *V2Add(Vector2 *a, Vector2 *b, Vector2 *c)
+{
+	(*c)[0] = (*a)[0]+(*b)[0];  (*c)[1] = (*a)[1]+(*b)[1];
+	return(c);
+} 
+
+/* normalizes the input vector and returns it */
+Vector2 *V2Normalize(Vector2 *v) 
+{
+double len = V2Length(v);
+	if (len != 0.0) { (*v)[0] /= len;  (*v)[1] /= len; }
+	return(v);
+}
+
+/* negates the input vector and returns it */
+Vector2 *V2Negate(Vector2 *v) 
+{
+	(*v)[0] = -(*v)[0];  (*v)[1] = -(*v)[1];
+	return(v);
+}
+
+ 
+/*
+ *  GenerateBezier :
+ *  Use least-squares method to find Bezier control points for region.
+ *
+ */
+static BezierCurve  GenerateBezier(Vector2 *d, int first, int last, double *uPrime, Vector2 tHat1, Vector2 tHat2)
+//    Vector2	*d;			/*  Array of digitized points	*/
+//    int		first, last;		/*  Indices defining region	*/
+//    double	*uPrime;		/*  Parameter values for region */
+//    Vector2	tHat1, tHat2;	/*  Unit tangents at endpoints	*/
+{
+    int 	i;
+    Vector2 	A[MAXPOINTS][2];	/* Precomputed rhs for eqn	*/
+    int 	nPts;			/* Number of pts in sub-curve */
+    double 	C[2][2];			/* Matrix C		*/
+    double 	X[2];			/* Matrix X			*/
+    double 	det_C0_C1,		/* Determinants of matrices	*/
+    	   	det_C0_X,
+	   		det_X_C1;
+    double 	alpha_l,		/* Alpha values, left and right	*/
+    	   	alpha_r;
+    Vector2 	tmp;			/* Utility variable		*/
+    BezierCurve	bezCurve;	/* RETURN bezier curve ctl pts	*/
+
+    bezCurve = (Vector2 *)malloc(4 * sizeof(Vector2));
+    nPts = last - first + 1;
+
+ 
+    /* Compute the A's	*/
+    for (i = 0; i < nPts; i++) {
+		Vector2		v1, v2;
+		v1 = tHat1;
+		v2 = tHat2;
+		V2Scale(&v1, B1(uPrime[i]));
+		V2Scale(&v2, B2(uPrime[i]));
+		A[i][0] = v1;
+		A[i][1] = v2;
+    }
+
+    /* Create the C and X matrices	*/
+    C[0][0] = 0.0;
+    C[0][1] = 0.0;
+    C[1][0] = 0.0;
+    C[1][1] = 0.0;
+    X[0]    = 0.0;
+    X[1]    = 0.0;
+
+    for (i = 0; i < nPts; i++) {
+        C[0][0] += V2Dot(&A[i][0], &A[i][0]);
+		C[0][1] += V2Dot(&A[i][0], &A[i][1]);
+/*					C[1][0] += V2Dot(&A[i][0], &A[i][1]);*/	
+		C[1][0] = C[0][1];
+		C[1][1] += V2Dot(&A[i][1], &A[i][1]);
+
+		tmp = V2SubII(d[first + i],
+	        V2AddII(
+	          V2ScaleIII(d[first], B0(uPrime[i])),
+		    	V2AddII(
+		      		V2ScaleIII(d[first], B1(uPrime[i])),
+		        			V2AddII(
+	                  		V2ScaleIII(d[last], B2(uPrime[i])),
+	                    		V2ScaleIII(d[last], B3(uPrime[i]))))));
+	
+
+	X[0] += V2Dot(&((A[i])[0]), &tmp);
+	X[1] += V2Dot(&((A[i])[1]), &tmp);
+    }
+
+    /* Compute the determinants of C and X	*/
+    det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
+    det_C0_X  = C[0][0] * X[1]    - C[0][1] * X[0];
+    det_X_C1  = X[0]    * C[1][1] - X[1]    * C[0][1];
+
+    /* Finally, derive alpha values	*/
+    if (det_C0_C1 == 0.0) {
+		det_C0_C1 = (C[0][0] * C[1][1]) * 10e-12;
+    }
+    alpha_l = det_X_C1 / det_C0_C1;
+    alpha_r = det_C0_X / det_C0_C1;
+
+
+    /*  If alpha negative, use the Wu/Barsky heuristic (see text) */
+	/* (if alpha is 0, you get coincident control points that lead to
+	 * divide by zero in any subsequent NewtonRaphsonRootFind() call. */
+    if (alpha_l < 1.0e-6 || alpha_r < 1.0e-6) {
+		double	dist = V2DistanceBetween2Points(&d[last], &d[first]) /
+					3.0;
+
+		bezCurve[0] = d[first];
+		bezCurve[3] = d[last];
+		V2Add(&(bezCurve[0]), V2Scale(&(tHat1), dist), &(bezCurve[1]));
+		V2Add(&(bezCurve[3]), V2Scale(&(tHat2), dist), &(bezCurve[2]));
+		return (bezCurve);
+    }
+
+    /*  First and last control points of the Bezier curve are */
+    /*  positioned exactly at the first and last data points */
+    /*  Control points 1 and 2 are positioned an alpha distance out */
+    /*  on the tangent vectors, left and right, respectively */
+    bezCurve[0] = d[first];
+    bezCurve[3] = d[last];
+    V2Add(&bezCurve[0], V2Scale(&tHat1, alpha_l), &bezCurve[1]);
+    V2Add(&bezCurve[3], V2Scale(&tHat2, alpha_r), &bezCurve[2]);
+    return (bezCurve);
+}
+
+
+/*
+ *  Reparameterize:
+ *	Given set of points and their parameterization, try to find
+ *   a better parameterization.
+ *
+ */
+static double *Reparameterize(Vector2 *d, int first, int last, double *u, BezierCurve bezCurve)
+//    Vector2	*d;			/*  Array of digitized points	*/
+//    int		first, last;		/*  Indices defining region	*/
+//    double	*u;			/*  Current parameter values	*/
+//    BezierCurve	bezCurve;	/*  Current fitted curve	*/
+{
+    int 	nPts = last-first+1;	
+    int 	i;
+    double	*uPrime;		/*  New parameter values	*/
+
+    uPrime = (double *)malloc(nPts * sizeof(double));
+    for (i = first; i <= last; i++) {
+		uPrime[i-first] = NewtonRaphsonRootFind(bezCurve, d[i], u[i-
+					first]);
+    }
+    return (uPrime);
+}
+
+
+
+/*
+ *  NewtonRaphsonRootFind :
+ *	Use Newton-Raphson iteration to find better root.
+ */
+static double NewtonRaphsonRootFind(BezierCurve Q, Vector2 P, double u)
+//    BezierCurve	Q;			/*  Current fitted curve	*/
+//    Vector2 		P;		/*  Digitized point		*/
+//    double 		u;		/*  Parameter value for "P"	*/
+{
+    double 		numerator, denominator;
+    Vector2 		Q1[3], Q2[2];	/*  Q' and Q''			*/
+    Vector2		Q_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q''	*/
+    double 		uPrime;		/*  Improved u			*/
+    int 		i;
+    
+    /* Compute Q(u)	*/
+    Q_u = BezierII(3, Q, u);
+    
+    /* Generate control vertices for Q'	*/
+    for (i = 0; i <= 2; i++) {
+		Q1[i][0] = (Q[i+1][0] - Q[i][0]) * 3.0;
+		Q1[i][1] = (Q[i+1][1] - Q[i][1]) * 3.0;
+    }
+    
+    /* Generate control vertices for Q'' */
+    for (i = 0; i <= 1; i++) {
+		Q2[i][0] = (Q1[i+1][0] - Q1[i][0]) * 2.0;
+		Q2[i][1] = (Q1[i+1][1] - Q1[i][1]) * 2.0;
+    }
+    
+    /* Compute Q'(u) and Q''(u)	*/
+    Q1_u = BezierII(2, Q1, u);
+    Q2_u = BezierII(1, Q2, u);
+    
+    /* Compute f(u)/f'(u) */
+    numerator = (Q_u[0] - P[0]) * (Q1_u[0]) + (Q_u[1] - P[1]) * (Q1_u[1]);
+    denominator = (Q1_u[0]) * (Q1_u[0]) + (Q1_u[1]) * (Q1_u[1]) +
+		      	  (Q_u[0] - P[0]) * (Q2_u[0]) + (Q_u[1] - P[1]) * (Q2_u[1]);
+    
+    /* u = u - f(u)/f'(u) */
+    if(denominator == 0) // FIXME
+      return u;
+    uPrime = u - (numerator/denominator);
+    return (uPrime);
+}
+
+	
+		       
+/*
+ *  Bezier :
+ *  	Evaluate a Bezier curve at a particular parameter value
+ * 
+ */
+static Vector2 BezierII(int degree, Vector2 *V, double t)
+//    int		degree;		/* The degree of the bezier curve	*/
+//    Vector2 	*V;		/* Array of control points		*/
+//    double 	t;		/* Parametric value to find point for	*/
+{
+    int 	i, j;		
+    Vector2 	Q;	        /* Point on curve at parameter t	*/
+    Vector2 	*Vtemp;		/* Local copy of control points		*/
+
+    /* Copy array	*/
+    Vtemp = (Vector2 *)malloc((unsigned)((degree+1) 
+				* sizeof (Vector2)));
+    for (i = 0; i <= degree; i++) {
+		Vtemp[i] = V[i];
+    }
+
+    /* Triangle computation	*/
+    for (i = 1; i <= degree; i++) {	
+		for (j = 0; j <= degree-i; j++) {
+	    	Vtemp[j][0] = (1.0 - t) * Vtemp[j][0] + t * Vtemp[j+1][0];
+	    	Vtemp[j][1] = (1.0 - t) * Vtemp[j][1] + t * Vtemp[j+1][1];
+		}
+    }
+
+    Q = Vtemp[0];
+    free((void *)Vtemp);
+    return Q;
+}
+
+
+/*
+ *  B0, B1, B2, B3 :
+ *	Bezier multipliers
+ */
+static double B0(double u)
+{
+    double tmp = 1.0 - u;
+    return (tmp * tmp * tmp);
+}
+
+
+static double B1(double u)
+{
+    double tmp = 1.0 - u;
+    return (3 * u * (tmp * tmp));
+}
+
+static double B2(double u)
+{
+    double tmp = 1.0 - u;
+    return (3 * u * u * tmp);
+}
+
+static double B3(double u)
+{
+    return (u * u * u);
+}
+
+
+
+/*
+ * ComputeLeftTangent, ComputeRightTangent, ComputeCenterTangent :
+ *Approximate unit tangents at endpoints and "center" of digitized curve
+ */
+static Vector2 ComputeLeftTangent(Vector2 *d, int end)
+//    Vector2	*d;			/*  Digitized points*/
+//    int		end;		/*  Index to "left" end of region */
+{
+    Vector2	tHat1;
+    tHat1 = V2SubII(d[end+1], d[end]);
+    tHat1 = *V2Normalize(&tHat1);
+    return tHat1;
+}
+
+static Vector2 ComputeRightTangent(Vector2 *d, int end)
+//    Vector2	*d;			/*  Digitized points		*/
+//    int		end;		/*  Index to "right" end of region */
+{
+    Vector2	tHat2;
+    tHat2 = V2SubII(d[end-1], d[end]);
+    tHat2 = *V2Normalize(&tHat2);
+    return tHat2;
+}
+
+static Vector2 ComputeCenterTangent(Vector2 *d, int center)
+//    Vector2	*d;			/*  Digitized points			*/
+//    int		center;		/*  Index to point inside region	*/
+{
+    Vector2	V1, V2, tHatCenter;
+
+    V1 = V2SubII(d[center-1], d[center]);
+    V2 = V2SubII(d[center], d[center+1]);
+    tHatCenter[0] = (V1[0] + V2[0])/2.0;
+    tHatCenter[1] = (V1[1] + V2[1])/2.0;
+    tHatCenter = *V2Normalize(&tHatCenter);
+    return tHatCenter;
+}
+
+
+/*
+ *  ChordLengthParameterize :
+ *	Assign parameter values to digitized points 
+ *	using relative distances between points.
+ */
+static double *ChordLengthParameterize(Vector2 *d, int first, int last)
+//    Vector2	*d;			/* Array of digitized points */
+//    int		first, last;		/*  Indices defining region	*/
+{
+    int		i;	
+    double	*u;			/*  Parameterization		*/
+
+    u = (double *)malloc((unsigned)(last-first+1) * sizeof(double));
+
+    u[0] = 0.0;
+    for (i = first+1; i <= last; i++) {
+		u[i-first] = u[i-first-1] +
+	  			V2DistanceBetween2Points(&d[i], &d[i-1]);
+    }
+
+    for (i = first + 1; i <= last; i++) {
+		u[i-first] = u[i-first] / u[last-first];
+    }
+
+    return(u);
+}
+
+
+
+
+/*
+ *  ComputeMaxError :
+ *	Find the maximum squared distance of digitized points
+ *	to fitted curve.
+*/
+static double ComputeMaxError(Vector2 *d, int first, int last, BezierCurve bezCurve, double *u, int *splitPoint)
+//    Vector2	*d;			/*  Array of digitized points	*/
+//    int		first, last;		/*  Indices defining region	*/
+//    BezierCurve	bezCurve;		/*  Fitted Bezier curve		*/
+//    double	*u;			/*  Parameterization of points	*/
+//    int		*splitPoint;		/*  Point of maximum error	*/
+{
+    int		i;
+    double	maxDist;		/*  Maximum error		*/
+    double	dist;		/*  Current error		*/
+    Vector2	P;			/*  Point on curve		*/
+    Vector2	v;			/*  Vector from point to curve	*/
+
+    *splitPoint = (last - first + 1)/2;
+    maxDist = 0.0;
+    for (i = first + 1; i < last; i++) {
+		P = BezierII(3, bezCurve, u[i-first]);
+		v = V2SubII(P, d[i]);
+		dist = V2SquaredLength(&v);
+		if (dist >= maxDist) {
+	    	maxDist = dist;
+	    	*splitPoint = i;
+		}
+    }
+    return (maxDist);
+}
+static Vector2 V2AddII(Vector2 a, Vector2 b)
+{
+    Vector2	c;
+    c[0] = a[0] + b[0];  c[1] = a[1] + b[1];
+    return (c);
+}
+static Vector2 V2ScaleIII(Vector2 v, double s)
+{
+    Vector2 result;
+    result[0] = v[0] * s; result[1] = v[1] * s;
+    return (result);
+}
+
+static Vector2 V2SubII(Vector2 a, Vector2 b)
+{
+    Vector2	c;
+    c[0] = a[0] - b[0]; c[1] = a[1] - b[1];
+    return (c);
+}
+
+#ifdef __cplusplus
+}
+#endif
+
+
+//------------------------- WRAPPER -----------------------------//
+
+FitCurveWrapper::FitCurveWrapper()
+{
+}
+
+FitCurveWrapper::~FitCurveWrapper()
+{
+  _vertices.clear();
+}
+
+void FitCurveWrapper::DrawBezierCurve(int n, Vector2 *curve )
+{
+	for(int i=0; i<n+1; ++i)
+    _vertices.push_back(curve[i]);
+}
+
+void FitCurveWrapper::FitCurve(vector<Vec2d>& data, vector<Vec2d>& oCurve, double error)
+{
+  int size = data.size();
+  Vector2 *d = new Vector2[size];
+  for(int i=0; i<size; ++i)
+  {
+    d[i][0] = data[i][0];
+    d[i][1] = data[i][1];
+  }
+
+  FitCurve(d,size,error);
+
+  // copy results
+  for(vector<Vector2>::iterator v=_vertices.begin(), vend=_vertices.end();
+  v!=vend;
+  ++v)
+  {
+    oCurve.push_back(Vec2d(v->x(), v->y())) ;
+  }
+  
+}
+
+void FitCurveWrapper::FitCurve(Vector2 *d, int nPts, double error)
+{
+    Vector2 tHat1, tHat2;	/*  Unit tangent vectors at endpoints */
+
+    tHat1 = ComputeLeftTangent(d, 0);
+    tHat2 = ComputeRightTangent(d, nPts - 1);
+    FitCubic(d, 0, nPts - 1, tHat1, tHat2, error);
+}
+
+void FitCurveWrapper::FitCubic(Vector2 *d, int first, int last, Vector2 tHat1, Vector2 tHat2, double error)
+{
+    BezierCurve	bezCurve; /*Control points of fitted Bezier curve*/
+    double	*u;		/*  Parameter values for point  */
+    double	*uPrime;	/*  Improved parameter values */
+    double	maxError;	/*  Maximum fitting error	 */
+    int		splitPoint;	/*  Point to split point set at	 */
+    int		nPts;		/*  Number of points in subset  */
+    double	iterationError; /*Error below which you try iterating  */
+    int		maxIterations = 4; /*  Max times to try iterating  */
+    Vector2	tHatCenter;   	/* Unit tangent vector at splitPoint */
+    int		i;		
+
+    iterationError = error * error;
+    nPts = last - first + 1;
+
+    /*  Use heuristic if region only has two points in it */
+    if (nPts == 2) {
+	    double dist = V2DistanceBetween2Points(&d[last], &d[first]) / 3.0;
+
+		bezCurve = (Vector2 *)malloc(4 * sizeof(Vector2));
+		bezCurve[0] = d[first];
+		bezCurve[3] = d[last];
+		V2Add(&bezCurve[0], V2Scale(&tHat1, dist), &bezCurve[1]);
+		V2Add(&bezCurve[3], V2Scale(&tHat2, dist), &bezCurve[2]);
+		DrawBezierCurve(3, bezCurve);
+		free((void *)bezCurve);
+		return;
+    }
+
+    /*  Parameterize points, and attempt to fit curve */
+    u = ChordLengthParameterize(d, first, last);
+    bezCurve = GenerateBezier(d, first, last, u, tHat1, tHat2);
+
+    /*  Find max deviation of points to fitted curve */
+    maxError = ComputeMaxError(d, first, last, bezCurve, u, &splitPoint);
+    if (maxError < error) {
+		DrawBezierCurve(3, bezCurve);
+		free((void *)u);
+		free((void *)bezCurve);
+		return;
+    }
+
+
+    /*  If error not too large, try some reparameterization  */
+    /*  and iteration */
+    if (maxError < iterationError) {
+		for (i = 0; i < maxIterations; i++) {
+	    	uPrime = Reparameterize(d, first, last, u, bezCurve);
+	    	bezCurve = GenerateBezier(d, first, last, uPrime, tHat1, tHat2);
+	    	maxError = ComputeMaxError(d, first, last,
+				       bezCurve, uPrime, &splitPoint);
+	    	if (maxError < error) {
+			DrawBezierCurve(3, bezCurve);
+			free((void *)u);
+			free((void *)bezCurve);
+			return;
+	    }
+	    free((void *)u);
+	    u = uPrime;
+	}
+    }
+
+    /* Fitting failed -- split at max error point and fit recursively */
+    free((void *)u);
+    free((void *)bezCurve);
+    tHatCenter = ComputeCenterTangent(d, splitPoint);
+    FitCubic(d, first, splitPoint, tHat1, tHatCenter, error);
+    V2Negate(&tHatCenter);
+    FitCubic(d, splitPoint, last, tHatCenter, tHat2, error);
+
+}
+
diff --git a/source/blender/freestyle/intern/geometry/FitCurve.h b/source/blender/freestyle/intern/geometry/FitCurve.h
new file mode 100755
index 00000000000..ed7cbe34780
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/FitCurve.h
@@ -0,0 +1,101 @@
+//
+//  Filename         : FitCurve.h
+//  Author(s)        : Stephane Grabli
+//  Purpose          : An Algorithm for Automatically Fitting Digitized Curves
+//                       by Philip J. Schneider
+//                       from "Graphics Gems", Academic Press, 1990
+//  Date of creation : 06/06/2003
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef  FITCURVE_H
+# define FITCURVE_H
+
+#include <vector>
+#include "../system/FreestyleConfig.h"
+#include "Geom.h"
+
+using namespace Geometry;
+
+typedef struct Point2Struct {	/* 2d point */
+	double coordinates[2];
+  Point2Struct() {coordinates[0]=0;coordinates[1]=0;}
+  inline double operator[](const int i) const 
+  {
+    return coordinates[i];
+  }
+  inline double& operator[](const int i)
+  {
+    return coordinates[i];
+  }
+  inline double x() const {return coordinates[0];}
+  inline double y() const {return coordinates[1];}
+	} Point2;
+
+typedef Point2 Vector2;
+
+
+
+class LIB_GEOMETRY_EXPORT FitCurveWrapper
+{
+private:
+  std::vector<Vector2> _vertices;
+
+public:
+  FitCurveWrapper();
+  ~FitCurveWrapper();
+
+  /*! Fits a set of 2D data points to a set of Bezier Curve segments
+   *    data
+   *      Input data points
+   *    oCurve
+   *      Control points of the sets of bezier curve segments.
+   *      Each segment is made of 4 points (polynomial degree of curve = 3)
+   *    error
+   *      max error tolerance between resulting curve and input data 
+   */
+  void FitCurve(std::vector<Vec2d>& data, std::vector<Vec2d>& oCurve, double error);
+protected:
+    /* Vec2d	*d;			  Array of digitized points	  */
+    /* int		nPts;		  Number of digitized points	*/
+    /* double	error;		User-defined error squared	*/
+  void FitCurve(Vector2 *d, int nPts, double error);
+    
+  /*! Draws a Bezier curve segment
+   *  n
+   *    degree of curve (=3)
+   *  curve
+   *    bezier segments control points
+   */
+  void DrawBezierCurve(int n, Vector2 *curve);
+  
+  /* Vec2d	*d;			  Array of digitized points */
+  /* int		first, last;	 Indices of first and last pts in region */
+  /* Vec2d tHat1, tHat2;	 Unit tangent vectors at endpoints */
+  /* double	error;		  User-defined error squared	   */
+  void FitCubic(Vector2 *d, int first, int last, Vector2 tHat1, Vector2 tHat2, double error);
+    
+};
+
+#endif // FITCURVE_H
diff --git a/source/blender/freestyle/intern/geometry/Geom.h b/source/blender/freestyle/intern/geometry/Geom.h
new file mode 100755
index 00000000000..ac94213fe98
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/Geom.h
@@ -0,0 +1,78 @@
+//
+//  Filename         : Geom.h
+//  Author(s)        : Sylvain Paris
+//                     Emmanuel Turquin
+//                     Stephane Grabli 
+//  Purpose          : Vectors and Matrices (useful type definitions)
+//  Date of creation : 20/05/2003
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef  GEOM_H
+# define GEOM_H
+
+# include "VecMat.h"
+# include "../system/Precision.h"
+
+namespace Geometry {
+
+  typedef VecMat::Vec2<unsigned>	Vec2u;
+  typedef VecMat::Vec2<int>		Vec2i;
+  typedef VecMat::Vec2<float>		Vec2f;
+  typedef VecMat::Vec2<double>		Vec2d;
+  typedef VecMat::Vec2<real>		Vec2r;
+
+  typedef VecMat::Vec3<unsigned>	Vec3u;
+  typedef VecMat::Vec3<int>		Vec3i;
+  typedef VecMat::Vec3<float>		Vec3f;
+  typedef VecMat::Vec3<double>		Vec3d;
+  typedef VecMat::Vec3<real>		Vec3r;
+
+  typedef VecMat::HVec3<unsigned>	HVec3u;
+  typedef VecMat::HVec3<int>		HVec3i;
+  typedef VecMat::HVec3<float>		HVec3f;
+  typedef VecMat::HVec3<double>		HVec3d;
+  typedef VecMat::HVec3<real>		HVec3r;
+
+  typedef VecMat::SquareMatrix<unsigned, 2>	Matrix22u;
+  typedef VecMat::SquareMatrix<int, 2>		Matrix22i;
+  typedef VecMat::SquareMatrix<float, 2>	Matrix22f;
+  typedef VecMat::SquareMatrix<double, 2>	Matrix22d;
+  typedef VecMat::SquareMatrix<real, 2>		Matrix22r;
+
+  typedef VecMat::SquareMatrix<unsigned, 3>	Matrix33u;
+  typedef VecMat::SquareMatrix<int, 3>		Matrix33i;
+  typedef VecMat::SquareMatrix<float, 3>	Matrix33f;
+  typedef VecMat::SquareMatrix<double, 3>	Matrix33d;
+  typedef VecMat::SquareMatrix<real, 3>		Matrix33r;
+
+  typedef VecMat::SquareMatrix<unsigned, 4>	Matrix44u;
+  typedef VecMat::SquareMatrix<int, 4>		Matrix44i;
+  typedef VecMat::SquareMatrix<float, 4>	Matrix44f;
+  typedef VecMat::SquareMatrix<double, 4>	Matrix44d;
+  typedef VecMat::SquareMatrix<real, 4>		Matrix44r;
+
+} // end of namespace Geometry
+
+#endif // GEOM_H
diff --git a/source/blender/freestyle/intern/geometry/GeomCleaner.cpp b/source/blender/freestyle/intern/geometry/GeomCleaner.cpp
new file mode 100755
index 00000000000..c148c521a46
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/GeomCleaner.cpp
@@ -0,0 +1,240 @@
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+//#if defined(__GNUC__) && (__GNUC__ >= 3)
+//// hash_map is not part of the C++ standard anymore;
+//// hash_map.h has been kept though for backward compatibility
+//# include <hash_map.h>
+//#else
+//# include <hash_map>
+//#endif
+
+#include <stdio.h>
+#include <list>
+#include <map>
+#include "../system/TimeUtils.h"
+#include "GeomCleaner.h"
+
+using namespace std;
+
+
+void GeomCleaner::SortIndexedVertexArray( const float *iVertices, unsigned iVSize, 
+                                          const unsigned *iIndices, unsigned iISize, 
+                                          real **oVertices, 
+                                          unsigned **oIndices)
+{
+  // First, we build a list of IndexVertex:
+  list<IndexedVertex> indexedVertices;
+  unsigned i;
+  for(i=0; i<iVSize; i+= 3)
+  {
+    indexedVertices.push_back(IndexedVertex(Vec3r(iVertices[i], iVertices[i+1], iVertices[i+2]), i/3));
+  }
+
+  // q-sort
+  indexedVertices.sort();
+
+  // build the indices mapping array:
+  unsigned *mapIndices = new unsigned[iVSize/3];
+  *oVertices = new real[iVSize];
+  list<IndexedVertex>::iterator iv;
+  unsigned newIndex = 0;
+  unsigned vIndex = 0;
+  for(iv=indexedVertices.begin(); iv!=indexedVertices.end(); iv++)
+  {
+    // Build the final results:
+    (*oVertices)[vIndex] = iv->x();
+    (*oVertices)[vIndex+1] = iv->y();
+    (*oVertices)[vIndex+2] = iv->z();
+
+    mapIndices[iv->index()] = newIndex;
+    newIndex++;
+    vIndex+=3;
+  }
+
+  
+  // Build the final index array:
+  *oIndices  = new unsigned[iISize];
+  for(i=0; i<iISize; i++)
+  {
+    (*oIndices)[i] = 3*mapIndices[iIndices[i]/3];
+  }
+
+  delete [] mapIndices;
+}
+
+void GeomCleaner::CompressIndexedVertexArray(const real *iVertices, unsigned iVSize, 
+                                             const unsigned *iIndices, unsigned iISize, 
+                                             real **oVertices, unsigned *oVSize, 
+                                             unsigned **oIndices)
+{
+  // First, we build a list of IndexVertex:
+  vector<Vec3r> vertices;
+  unsigned i;
+  for(i=0; i<iVSize; i+= 3)
+  {
+    vertices.push_back(Vec3r(iVertices[i], iVertices[i+1], iVertices[i+2]));
+  }
+
+  unsigned *mapVertex = new unsigned[iVSize];
+  vector<Vec3r>::iterator v = vertices.begin();
+  
+  vector<Vec3r> compressedVertices;
+  Vec3r previous = *v;
+  mapVertex[0] = 0;
+  compressedVertices.push_back(vertices.front());
+
+  v++;
+  Vec3r current;
+  i=1;
+  for(; v!=vertices.end(); v++)
+  {
+    current = *v;
+    if(current == previous)
+      mapVertex[i] = compressedVertices.size()-1;
+    else
+    {
+      compressedVertices.push_back(current);
+      mapVertex[i] = compressedVertices.size()-1;
+    }  
+    previous = current;
+    i++;
+  }
+
+  // Builds the resulting vertex array:
+  *oVSize = 3*compressedVertices.size();
+  *oVertices = new real [*oVSize];
+  i=0;
+  for(v=compressedVertices.begin(); v!=compressedVertices.end(); v++)
+  {
+    (*oVertices)[i] = (*v)[0];
+    (*oVertices)[i+1] = (*v)[1];
+    (*oVertices)[i+2] = (*v)[2];
+    i += 3;
+  }
+
+  // Map the index array:
+  *oIndices = new unsigned[iISize];
+  for(i=0; i<iISize; i++)
+  {
+    (*oIndices)[i] = 3*mapVertex[iIndices[i]/3];
+  }
+
+  delete [] mapVertex;
+}
+
+void GeomCleaner::SortAndCompressIndexedVertexArray(const float *iVertices, unsigned iVSize, 
+                                                    const unsigned *iIndices, unsigned iISize, 
+                                                    real **oVertices, unsigned *oVSize, 
+                                                    unsigned **oIndices)
+{
+
+  // tmp arrays used to store the sorted data:
+  real *tmpVertices;
+  unsigned   *tmpIndices;
+
+  Chronometer chrono;
+  // Sort data
+  chrono.start();
+  GeomCleaner::SortIndexedVertexArray(iVertices, iVSize, 
+                                      iIndices, iISize, 
+                                      &tmpVertices, &tmpIndices
+                                      );
+  printf("Sorting: %lf\n", chrono.stop());
+
+  // compress data
+  chrono.start();
+  GeomCleaner::CompressIndexedVertexArray(tmpVertices, iVSize, 
+                                          tmpIndices, iISize, 
+                                          oVertices, oVSize,
+                                          oIndices);
+  printf("Merging: %lf\n", chrono.stop());
+
+  // deallocates memory:
+  delete [] tmpVertices;
+  delete [] tmpIndices;
+}
+
+/*! Defines a hash table used for searching the Cells */
+struct GeomCleanerHasher{
+#define _MUL 950706376UL
+#define _MOD 2147483647UL
+  inline size_t operator() (const Vec3r& p) const {
+    size_t res = ((unsigned long) (p[0] * _MUL)) % _MOD;
+    res = ((res + (unsigned long) (p[1]) * _MUL)) % _MOD;
+    return ((res +(unsigned long) (p[2]) * _MUL)) % _MOD;
+  }
+};
+
+void GeomCleaner::CleanIndexedVertexArray(const float *iVertices, unsigned iVSize, 
+					  const unsigned *iIndices, unsigned iISize, 
+					  real **oVertices, unsigned *oVSize, 
+					  unsigned **oIndices)
+{
+  typedef map<Vec3r, unsigned> cleanHashTable;
+  vector<Vec3r> vertices;
+  unsigned i;
+  for(i=0; i<iVSize; i+= 3)
+    vertices.push_back(Vec3r(iVertices[i], iVertices[i+1], iVertices[i+2]));
+  
+  cleanHashTable ht;
+  vector<unsigned> newIndices;
+  vector<Vec3r> newVertices;
+  
+  // elimination of needless points
+  unsigned currentIndex = 0;
+  vector<Vec3r>::const_iterator v = vertices.begin();
+  vector<Vec3r>::const_iterator end = vertices.end();
+  cleanHashTable::const_iterator found;
+  for(; v!=end; v++)
+  {
+    found = ht.find(*v);
+    if(found != ht.end())
+    {
+      // The vertex is already in the new array.
+      newIndices.push_back((*found).second);
+    }
+    else
+    {
+      newVertices.push_back(*v);
+      newIndices.push_back(currentIndex);
+      ht[*v] = currentIndex;
+      currentIndex++;
+    }
+  }
+
+  // creation of oVertices array:
+  *oVSize = 3*newVertices.size();
+  *oVertices = new real[*oVSize];
+  currentIndex = 0;
+  end = newVertices.end();
+  for(v=newVertices.begin(); v!=end ; v++)
+  {
+    (*oVertices)[currentIndex++] = (*v)[0];
+    (*oVertices)[currentIndex++] = (*v)[1];
+    (*oVertices)[currentIndex++] = (*v)[2];
+  }
+
+  // map new indices:
+  *oIndices = new unsigned[iISize];
+  for(i=0; i<iISize; i++)
+    (*oIndices)[i] = 3*newIndices[iIndices[i]/3];
+}
diff --git a/source/blender/freestyle/intern/geometry/GeomCleaner.h b/source/blender/freestyle/intern/geometry/GeomCleaner.h
new file mode 100755
index 00000000000..d78d90ccb4a
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/GeomCleaner.h
@@ -0,0 +1,219 @@
+//
+//  Filename         : GeomCleaner.h
+//  Author           : Stephane Grabli
+//  Purpose          : Class to define a cleaner of geometry providing 
+//                     a set of useful tools
+//  Date of creation : 04/03/2002
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef  GEOMCLEANER_H
+# define GEOMCLEANER_H
+
+# include "../system/FreestyleConfig.h"
+# include "Geom.h"
+
+using namespace Geometry;
+
+class LIB_GEOMETRY_EXPORT GeomCleaner
+{
+public:
+
+  inline GeomCleaner() {}
+  inline ~GeomCleaner() {}
+
+  /*! Sorts an array of Indexed vertices 
+   *    iVertices
+   *      Array of vertices to sort. It is organized as a
+   *      float series of vertex coordinates: XYZXYZXYZ...
+   *    iVSize
+   *      The size of iVertices array.
+   *    iIndices
+   *      The array containing the vertex indices (used to refer
+   *      to the vertex coordinates in an indexed face). Each 
+   *      element is an unsignedeger multiple of 3.
+   *    iISize
+   *      The size of iIndices array
+   *    oVertices
+   *      Output of sorted vertices. A vertex v1 precedes another one 
+   *      v2 in this array if v1.x<v2.x, or v1.x=v2.x && v1.y < v2.y 
+   *      or v1.x=v2.y && v1.y=v2.y && v1.z < v2.z.
+   *      The array is organized as a 3-float serie giving 
+   *      the vertices coordinates: XYZXYZXYZ...
+   *    oIndices
+   *      Output corresponding to the iIndices array but 
+   *      reorganized in order to match the sorted vertex array.
+   */    
+
+  static void SortIndexedVertexArray(const float *iVertices, unsigned iVSize, 
+                                     const unsigned *iIndices, unsigned iISize, 
+                                     real **oVertices, 
+                                     unsigned **oIndices);
+
+  /*! Compress a SORTED indexed vertex array by eliminating 
+   *  multiple appearing occurences of a single vertex.
+   *    iVertices
+   *      The SORTED vertex array to compress. It is organized as a
+   *      float series of vertex coordinates: XYZXYZXYZ...
+   *    iVSize
+   *      The size of iVertices array.
+   *    iIndices
+   *      The array containing the vertex indices (used to refer
+   *      to the vertex coordinates in an indexed face). Each 
+   *      element is an unsignedeger multiple of 3.
+   *    iISize
+   *      The size of iIndices array
+   *    oVertices
+   *      The vertex array, result of the compression.
+   *      The array is organized as a 3-float serie giving 
+   *      the vertices coordinates: XYZXYZXYZ...
+   *    oVSize
+   *      The size of oVertices.
+   *    oIndices
+   *      The indices array, reorganized to match the compressed 
+   *      oVertices array.
+   */
+
+  static void CompressIndexedVertexArray(const real *iVertices, unsigned iVSize, 
+                                         const unsigned *iIndices, unsigned iISize, 
+                                         real **oVertices, unsigned *oVSize, 
+                                         unsigned **oIndices);
+
+  /*! Sorts and compress an array of indexed vertices.
+   *    iVertices
+   *      The vertex array to sort then compress. It is organized as a
+   *      float series of vertex coordinates: XYZXYZXYZ...
+   *    iVSize
+   *      The size of iVertices array.
+   *    iIndices
+   *      The array containing the vertex indices (used to refer
+   *      to the vertex coordinates in an indexed face). Each 
+   *      element is an unsignedeger multiple of 3.
+   *    iISize
+   *      The size of iIndices array
+   *    oVertices
+   *      The vertex array, result of the sorting-compression.
+   *      The array is organized as a 3-float serie giving 
+   *      the vertices coordinates: XYZXYZXYZ...
+   *    oVSize
+   *      The size of oVertices.
+   *    oIndices
+   *      The indices array, reorganized to match the sorted and compressed 
+   *      oVertices array.
+   */
+
+  static void SortAndCompressIndexedVertexArray(const float *iVertices, unsigned iVSize, 
+                                                const unsigned *iIndices, unsigned iISize, 
+                                                real **oVertices, unsigned *oVSize, 
+                                                unsigned **oIndices);
+
+  /*! Cleans an indexed vertex array. (Identical to 
+   *  SortAndCompress except that we use here a hash
+   *  table to create the new array.)
+   *    iVertices
+   *      The vertex array to sort then compress. It is organized as a
+   *      float series of vertex coordinates: XYZXYZXYZ...
+   *    iVSize
+   *      The size of iVertices array.
+   *    iIndices
+   *      The array containing the vertex indices (used to refer
+   *      to the vertex coordinates in an indexed face). Each 
+   *      element is an unsignedeger multiple of 3.
+   *    iISize
+   *      The size of iIndices array
+   *    oVertices
+   *      The vertex array, result of the sorting-compression.
+   *      The array is organized as a 3-float serie giving 
+   *      the vertices coordinates: XYZXYZXYZ...
+   *    oVSize
+   *      The size of oVertices.
+   *    oIndices
+   *      The indices array, reorganized to match the sorted and compressed 
+   *      oVertices array.
+   */
+
+  static void CleanIndexedVertexArray(const float *iVertices, unsigned iVSize, 
+                                      const unsigned *iIndices, unsigned iISize, 
+                                      real **oVertices, unsigned *oVSize, 
+                                      unsigned **oIndices);
+};
+
+
+/*! Binary operators */
+//inline bool operator<(const IndexedVertex& iv1, const IndexedVertex& iv2);
+
+/*! Class Indexed Vertex. Used to represent 
+ *  an indexed vertex by storing the vertex 
+ *  coordinates as well as its index
+ */
+class IndexedVertex
+{
+public:
+
+private:
+  Vec3r      _Vector;
+  unsigned        _index;
+public:
+  inline IndexedVertex() {}
+  inline IndexedVertex(Vec3r iVector, unsigned iIndex)
+  {
+    _Vector = iVector;
+    _index = iIndex;
+  }
+  /*! accessors */
+  inline const Vec3r& vector() const {return _Vector;}
+  inline unsigned index() {return _index;}
+  inline real x() {return _Vector[0];}
+  inline real y() {return _Vector[1];}
+  inline real z() {return _Vector[2];}
+  
+  /*! modifiers */
+  inline void SetVector(const Vec3r& iVector) {_Vector = iVector;}
+  inline void SetIndex(unsigned iIndex) {_index = iIndex;}
+
+  /*! operators */
+  IndexedVertex& operator=(const IndexedVertex& iv)
+  {
+    _Vector = iv._Vector;
+    _index  = iv._index;
+    return *this;
+  }
+  inline real operator[](const unsigned i) {return _Vector[i];}
+  //friend inline bool operator<(const IndexedVertex& iv1, const IndexedVertex& iv2);
+  inline bool operator<(const IndexedVertex& v) const
+  {
+    return (_Vector < v._Vector);
+  }
+  inline bool operator==(const IndexedVertex& v)
+  {
+    return (_Vector == v._Vector);
+  }
+};
+
+//bool operator<(const IndexedVertex& iv1, const IndexedVertex& iv2)
+//{
+//  return iv1.operator<(iv2);
+//}
+
+#endif // GEOMCLEANER_H
diff --git a/source/blender/freestyle/intern/geometry/GeomUtils.cpp b/source/blender/freestyle/intern/geometry/GeomUtils.cpp
new file mode 100755
index 00000000000..fd36e81ca77
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/GeomUtils.cpp
@@ -0,0 +1,742 @@
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#include "GeomUtils.h"
+
+namespace GeomUtils {
+
+  // This internal procedure is defined below.
+  bool intersect2dSegPoly(Vec2r* seg,
+			  Vec2r* poly,
+			  unsigned n);
+
+  bool intersect2dSeg2dArea(const Vec2r& min,
+			    const Vec2r& max,
+			    const Vec2r& A,
+			    const Vec2r& B) {
+    Vec2r seg[2];
+    seg[0] = A;
+    seg[1] = B;
+
+    Vec2r poly[5];
+    poly[0][0] = min[0];
+    poly[0][1] = min[1];
+    poly[1][0] = max[0];
+    poly[1][1] = min[1];
+    poly[2][0] = max[0];
+    poly[2][1] = max[1];
+    poly[3][0] = min[0];
+    poly[3][1] = max[1];
+    poly[4][0] = min[0];
+    poly[4][1] = min[1];
+
+    return intersect2dSegPoly(seg, poly, 4);
+  }
+
+  bool include2dSeg2dArea(const Vec2r& min,
+			  const Vec2r& max,
+			  const Vec2r& A,
+			  const Vec2r& B) {
+    if((((max[0] > A[0])&&(A[0] > min[0]))&&((max[0] > B[0])&&(B[0] > min[0])))
+       && (((max[1] > A[1])&&(A[1] > min[1]))&&((max[1] > B[1])&&(B[1] > min[1]))))
+      return true;
+    return false;
+  }
+
+  intersection_test intersect2dSeg2dSeg(const Vec2r& p1,
+					const Vec2r& p2,
+					const Vec2r& p3,
+					const Vec2r& p4,
+					Vec2r& res) {
+    real a1, a2, b1, b2, c1, c2; // Coefficients of line eqns
+    real r1, r2, r3, r4;         // 'Sign' values
+    real denom, num;             // Intermediate values
+
+    // Compute a1, b1, c1, where line joining points p1 and p2
+    // is "a1 x  +  b1 y  +  c1  =  0".
+    a1 = p2[1] - p1[1];
+    b1 = p1[0] - p2[0];
+    c1 = p2[0] * p1[1] - p1[0] * p2[1];
+
+    // Compute r3 and r4.
+    r3 = a1 * p3[0] + b1 * p3[1] + c1;
+    r4 = a1 * p4[0] + b1 * p4[1] + c1;
+
+    // Check signs of r3 and r4.  If both point 3 and point 4 lie on
+    // same side of line 1, the line segments do not intersect.
+    if ( r3 != 0 && r4 != 0 && r3 * r4 > 0.0)
+      return (DONT_INTERSECT);
+
+    // Compute a2, b2, c2
+    a2 = p4[1] - p3[1];
+    b2 = p3[0] - p4[0];
+    c2 = p4[0] * p3[1] - p3[0] * p4[1];
+
+    // Compute r1 and r2
+    r1 = a2 * p1[0] + b2 * p1[1] + c2;
+    r2 = a2 * p2[0] + b2 * p2[1] + c2;
+
+    // Check signs of r1 and r2.  If both point 1 and point 2 lie
+    // on same side of second line segment, the line segments do
+    // not intersect.
+    if ( r1 != 0 && r2 != 0 && r1 * r2 > 0.0)
+      return (DONT_INTERSECT);
+
+    // Line segments intersect: compute intersection point.
+    denom = a1 * b2 - a2 * b1;
+    if (fabs(denom) < M_EPSILON)
+      return (COLINEAR);
+    
+    num = b1 * c2 - b2 * c1;
+    res[0] = num / denom;
+
+    num = a2 * c1 - a1 * c2;
+    res[1] = num / denom;
+
+    return (DO_INTERSECT);
+  }
+
+  intersection_test intersect2dLine2dLine(const Vec2r& p1,
+					  const Vec2r& p2,
+					  const Vec2r& p3,
+					  const Vec2r& p4,
+					  Vec2r& res) {
+    real a1, a2, b1, b2, c1, c2; // Coefficients of line eqns
+    real denom, num;             // Intermediate values
+
+    // Compute a1, b1, c1, where line joining points p1 and p2
+    // is "a1 x  +  b1 y  +  c1  =  0".
+    a1 = p2[1] - p1[1];
+    b1 = p1[0] - p2[0];
+    c1 = p2[0] * p1[1] - p1[0] * p2[1];
+
+    // Compute a2, b2, c2
+    a2 = p4[1] - p3[1];
+    b2 = p3[0] - p4[0];
+    c2 = p4[0] * p3[1] - p3[0] * p4[1];
+
+    // Line segments intersect: compute intersection point.
+    denom = a1 * b2 - a2 * b1;
+    if (fabs(denom) < M_EPSILON)
+      return (COLINEAR);
+
+    num = b1 * c2 - b2 * c1;
+    res[0] = num / denom;
+
+    num = a2 * c1 - a1 * c2;
+    res[1] = num / denom;
+
+    return (DO_INTERSECT);
+  }
+
+  intersection_test intersect2dSeg2dSegParametric(const Vec2r& p1,
+						  const Vec2r& p2,
+						  const Vec2r& p3,
+						  const Vec2r& p4,
+						  real& t,
+						  real& u) {
+    real a1, a2, b1, b2, c1, c2; // Coefficients of line eqns
+    real r1, r2, r3, r4;         // 'Sign' values
+    real denom, num;             // Intermediate values
+
+    // Compute a1, b1, c1, where line joining points p1 and p2
+    // is "a1 x  +  b1 y  +  c1  =  0".
+    a1 = p2[1] - p1[1];
+    b1 = p1[0] - p2[0];
+    c1 = p2[0] * p1[1] - p1[0] * p2[1];
+
+    // Compute r3 and r4.
+    r3 = a1 * p3[0] + b1 * p3[1] + c1;
+    r4 = a1 * p4[0] + b1 * p4[1] + c1;
+
+    // Check signs of r3 and r4.  If both point 3 and point 4 lie on
+    // same side of line 1, the line segments do not intersect.
+    if ( r3 != 0 && r4 != 0 && r3 * r4 > 0.0)
+      return (DONT_INTERSECT);
+
+    // Compute a2, b2, c2
+    a2 = p4[1] - p3[1];
+    b2 = p3[0] - p4[0];
+    c2 = p4[0] * p3[1] - p3[0] * p4[1];
+
+    // Compute r1 and r2
+    r1 = a2 * p1[0] + b2 * p1[1] + c2;
+    r2 = a2 * p2[0] + b2 * p2[1] + c2;
+
+    // Check signs of r1 and r2.  If both point 1 and point 2 lie
+    // on same side of second line segment, the line segments do
+    // not intersect.
+    if ( r1 != 0 && r2 != 0 && r1 * r2 > 0.0)
+      return (DONT_INTERSECT);
+
+    // Line segments intersect: compute intersection point.
+    denom = a1 * b2 - a2 * b1;
+    if (fabs(denom) < M_EPSILON)
+      return (COLINEAR);
+    
+    real d1, d2, e1;
+
+    d1 = p1[1] - p3[1];
+    d2 = p2[1] - p1[1];
+    e1 = p1[0] - p3[0];
+
+    num = -b2 * d1 - a2 * e1;
+    t = num / denom;
+
+    num = -b1 * d1 - a1 * e1;
+    u = num / denom;
+
+    return (DO_INTERSECT);
+  }
+
+  // AABB-triangle overlap test code
+  // by Tomas Akenine-Möller
+  // Function: int triBoxOverlap(real boxcenter[3],
+  //          real boxhalfsize[3],real triverts[3][3]);
+  // History:
+  //   2001-03-05: released the code in its first version
+  //   2001-06-18: changed the order of the tests, faster
+  //
+  // Acknowledgement: Many thanks to Pierre Terdiman for
+  // suggestions and discussions on how to optimize code.
+  // Thanks to David Hunt for finding a ">="-bug!
+
+#define X 0
+#define Y 1
+#define Z 2
+
+#define FINDMINMAX(x0, x1, x2, min, max) \
+  min = max = x0;    \
+  if(x1<min) min=x1; \
+  if(x1>max) max=x1; \
+  if(x2<min) min=x2; \
+  if(x2>max) max=x2;
+
+  //======================== X-tests ========================//
+#define AXISTEST_X01(a, b, fa, fb)                         \
+        p0 = a*v0[Y] - b*v0[Z];                            \
+        p2 = a*v2[Y] - b*v2[Z];                            \
+        if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
+        rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z];   \
+        if(min>rad || max<-rad) return 0;
+
+#define AXISTEST_X2(a, b, fa, fb)                          \
+        p0 = a*v0[Y] - b*v0[Z];                            \
+        p1 = a*v1[Y] - b*v1[Z];                            \
+        if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
+        rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z];   \
+        if(min>rad || max<-rad) return 0;
+
+  //======================== Y-tests ========================//
+#define AXISTEST_Y02(a, b, fa, fb)                         \
+        p0 = -a*v0[X] + b*v0[Z];                           \
+        p2 = -a*v2[X] + b*v2[Z];                           \
+        if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
+        rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z];   \
+        if(min>rad || max<-rad) return 0;
+
+#define AXISTEST_Y1(a, b, fa, fb)                          \
+        p0 = -a*v0[X] + b*v0[Z];                           \
+        p1 = -a*v1[X] + b*v1[Z];                           \
+        if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
+        rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z];   \
+        if(min>rad || max<-rad) return 0;
+
+  //======================== Z-tests ========================//
+#define AXISTEST_Z12(a, b, fa, fb)                         \
+        p1 = a*v1[X] - b*v1[Y];                            \
+        p2 = a*v2[X] - b*v2[Y];                            \
+        if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;} \
+        rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y];   \
+        if(min>rad || max<-rad) return 0;
+
+#define AXISTEST_Z0(a, b, fa, fb)                          \
+        p0 = a*v0[X] - b*v0[Y];                            \
+        p1 = a*v1[X] - b*v1[Y];                            \
+        if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
+        rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y];   \
+        if(min>rad || max<-rad) return 0;
+ 
+  // This internal procedure is defined below.
+  bool overlapPlaneBox(Vec3r& normal, real d, Vec3r& maxbox);
+
+  // Use separating axis theorem to test overlap between triangle and box
+  // need to test for overlap in these directions:
+  // 1) the {x,y,z}-directions (actually, since we use the AABB of the triangle
+  //    we do not even need to test these)
+  // 2) normal of the triangle
+  // 3) crossproduct(edge from tri, {x,y,z}-directin)
+  //    this gives 3x3=9 more tests
+  bool overlapTriangleBox(Vec3r& boxcenter,
+			  Vec3r& boxhalfsize,
+			  Vec3r triverts[3]) {
+    Vec3r v0, v1, v2, normal, e0, e1, e2;
+    real min, max, d, p0, p1, p2, rad, fex, fey, fez;  
+
+    // This is the fastest branch on Sun
+    // move everything so that the boxcenter is in (0, 0, 0)
+    v0 = triverts[0] - boxcenter;
+    v1 = triverts[1] - boxcenter;
+    v2 = triverts[2] - boxcenter;
+
+    // compute triangle edges
+    e0 = v1 - v0;
+    e1 = v2 - v1;
+    e2 = v0 - v2;
+
+    // Bullet 3:
+    // Do the 9 tests first (this was faster)
+    fex = fabs(e0[X]);
+    fey = fabs(e0[Y]);
+    fez = fabs(e0[Z]);
+    AXISTEST_X01(e0[Z], e0[Y], fez, fey);
+    AXISTEST_Y02(e0[Z], e0[X], fez, fex);
+    AXISTEST_Z12(e0[Y], e0[X], fey, fex);
+
+    fex = fabs(e1[X]);
+    fey = fabs(e1[Y]);
+    fez = fabs(e1[Z]);
+    AXISTEST_X01(e1[Z], e1[Y], fez, fey);
+    AXISTEST_Y02(e1[Z], e1[X], fez, fex);
+    AXISTEST_Z0(e1[Y], e1[X], fey, fex);
+  
+    fex = fabs(e2[X]);
+    fey = fabs(e2[Y]);
+    fez = fabs(e2[Z]);
+    AXISTEST_X2(e2[Z], e2[Y], fez, fey);
+    AXISTEST_Y1(e2[Z], e2[X], fez, fex);
+    AXISTEST_Z12(e2[Y], e2[X], fey, fex);
+
+    // Bullet 1:
+    // first test overlap in the {x,y,z}-directions
+    // find min, max of the triangle each direction, and test for overlap in
+    // that direction -- this is equivalent to testing a minimal AABB around
+    // the triangle against the AABB
+
+    // test in X-direction
+    FINDMINMAX(v0[X], v1[X], v2[X], min, max);
+    if (min > boxhalfsize[X] || max < -boxhalfsize[X])
+      return false;
+
+    // test in Y-direction
+    FINDMINMAX(v0[Y], v1[Y], v2[Y], min, max);
+    if (min > boxhalfsize[Y] || max < -boxhalfsize[Y])
+      return false;
+
+    // test in Z-direction
+    FINDMINMAX(v0[Z], v1[Z], v2[Z], min, max);
+    if(min > boxhalfsize[Z] || max < -boxhalfsize[Z])
+      return false;
+
+    // Bullet 2:
+    // test if the box intersects the plane of the triangle
+    // compute plane equation of triangle: normal * x + d = 0
+    normal = e0 ^ e1;
+    d = -(normal * v0); // plane eq: normal.x + d = 0
+    if (!overlapPlaneBox(normal, d, boxhalfsize))
+      return false;
+  
+    return true; // box and triangle overlaps
+  }
+
+  // Fast, Minimum Storage Ray-Triangle Intersection
+  //
+  // Tomas Möller
+  // Prosolvia Clarus AB
+  // Sweden
+  // tompa@clarus.se
+  //
+  // Ben Trumbore
+  // Cornell University
+  // Ithaca, New York
+  // wbt@graphics.cornell.edu
+
+  bool intersectRayTriangle(Vec3r& orig, Vec3r& dir,
+			    Vec3r& v0, Vec3r& v1, Vec3r& v2,
+			    real& t, real& u, real& v, real epsilon) {
+    Vec3r edge1, edge2, tvec, pvec, qvec;
+    real det, inv_det;
+
+    // find vectors for two edges sharing v0
+    edge1 = v1 - v0;
+    edge2 = v2 - v0;
+
+    // begin calculating determinant - also used to calculate U parameter
+    pvec = dir ^ edge2;
+
+    // if determinant is near zero, ray lies in plane of triangle
+    det = edge1 * pvec;
+
+    // calculate distance from v0 to ray origin
+    tvec = orig - v0;
+    inv_det = 1.0 / det;
+   
+    qvec = tvec ^ edge1;
+      
+    if (det > epsilon) {
+      u = tvec * pvec;
+      if (u < 0.0 || u > det)
+	return false;
+            
+      // calculate V parameter and test bounds
+      v = dir * qvec;
+      if (v < 0.0 || u + v > det)
+	return false;
+    }
+    else if(det < -epsilon) {
+      // calculate U parameter and test bounds
+      u = tvec * pvec;
+      if (u > 0.0 || u < det)
+	return false;
+      
+      // calculate V parameter and test bounds
+      v = dir * qvec;
+      if (v > 0.0 || u + v < det)
+	return false;
+    }
+    else
+      return false;  // ray is parallell to the plane of the triangle
+
+    u *= inv_det;
+    v *= inv_det;
+    t = (edge2 * qvec) * inv_det;
+
+    return true;
+  }
+
+  // Intersection between plane and ray, adapted from Graphics Gems, Didier Badouel
+  intersection_test intersectRayPlane(Vec3r& orig, Vec3r& dir,
+				      Vec3r& norm, real d,
+				      real& t,
+				      real epsilon) { 
+    real denom = norm * dir;
+
+    if(fabs(denom) <= epsilon) { // plane and ray are parallel
+      if(fabs((norm * orig) + d) <= epsilon)
+        return COINCIDENT; // plane and ray are coincident
+      else
+        return COLINEAR;
+    }
+
+    t = -(d + (norm * orig)) / denom;
+
+    if (t < 0.0f)
+      return DONT_INTERSECT;
+
+    return DO_INTERSECT;
+  }
+
+    bool intersectRayBBox(const Vec3r& orig, const Vec3r& dir,   // ray origin and direction
+        const Vec3r& boxMin, const Vec3r& boxMax, // the bbox
+        real t0, real t1,
+        real& tmin, real& tmax,                  // I0=orig+tmin*dir is the first intersection, I1=orig+tmax*dir is the second intersection
+        real epsilon){
+
+            float tymin, tymax, tzmin, tzmax;
+            Vec3r inv_direction(1.0/dir[0], 1.0/dir[1], 1.0/dir[2]);
+            int sign[3];
+            sign[0] = (inv_direction.x() < 0);
+            sign[1] = (inv_direction.y() < 0);
+            sign[2] = (inv_direction.z() < 0);
+
+            Vec3r bounds[2];
+            bounds[0] = boxMin;
+            bounds[1] = boxMax;
+
+            tmin = (bounds[sign[0]].x() - orig.x()) * inv_direction.x();
+            tmax = (bounds[1-sign[0]].x() - orig.x()) * inv_direction.x();
+            tymin = (bounds[sign[1]].y() - orig.y()) * inv_direction.y();
+            tymax = (bounds[1-sign[1]].y() - orig.y()) * inv_direction.y();
+            if ( (tmin > tymax) || (tymin > tmax) )
+                return false;
+            if (tymin > tmin)
+                tmin = tymin;
+            if (tymax < tmax)
+                tmax = tymax;
+            tzmin = (bounds[sign[2]].z() - orig.z()) * inv_direction.z();
+            tzmax = (bounds[1-sign[2]].z() - orig.z()) * inv_direction.z();
+            if ( (tmin > tzmax) || (tzmin > tmax) )
+                return false;
+            if (tzmin > tmin)
+                tmin = tzmin;
+            if (tzmax < tmax)
+                tmax = tzmax;
+            return ( (tmin < t1) && (tmax > t0) );
+        }
+
+  // Checks whether 3D points p lies inside or outside of the triangle ABC
+  bool includePointTriangle(Vec3r& P,
+			    Vec3r& A,
+			    Vec3r& B,
+			    Vec3r& C) {
+    Vec3r AB(B - A);
+    Vec3r BC(C - B);
+    Vec3r CA(A - C);
+    Vec3r AP(P - A);
+    Vec3r BP(P - B);
+    Vec3r CP(P - C);
+
+    Vec3r N(AB ^ BC); // triangle's normal
+
+    N.normalize();
+
+    Vec3r J(AB ^ AP), K(BC ^ BP), L(CA ^ CP);
+    J.normalize();
+    K.normalize();
+    L.normalize();
+
+    if(J * N < 0)
+      return false; // on the right of AB
+
+    if(K * N < 0)
+      return false; // on the right of BC
+
+    if(L * N < 0)
+      return false; // on the right of CA
+
+    return true;
+  }
+
+  void transformVertex(const Vec3r& vert,
+		       const Matrix44r& matrix,
+		       Vec3r& res) {
+    HVec3r hvert(vert), res_tmp;    
+    real scale;
+    for (unsigned j = 0; j < 4; j++) {
+      scale = hvert[j];
+      for (unsigned i = 0; i < 4; i++)
+        res_tmp[i] += matrix(i, j) * scale;
+    }
+    
+    res[0] = res_tmp.x();
+    res[1] = res_tmp.y();
+    res[2] = res_tmp.z();
+  }
+  
+  void transformVertices(const vector<Vec3r>& vertices,
+			 const Matrix44r& trans,
+			 vector<Vec3r>& res) {
+    for (vector<Vec3r>::const_iterator v = vertices.begin();
+	 v != vertices.end();
+	 v++) {
+      Vec3r *res_tmp = new Vec3r;
+      transformVertex(*v, trans, *res_tmp);
+      res.push_back(*res_tmp);
+    }
+  }
+
+  Vec3r rotateVector(const Matrix44r& mat, const Vec3r& v) {
+    Vec3r res;
+    for (unsigned i = 0; i < 3; i++) {
+      res[i] = 0;
+      for (unsigned j = 0; j < 3; j++)
+	res[i] += mat(i, j) * v[j];
+    }
+    res.normalize();
+    return res;
+  }
+
+  // This internal procedure is defined below.
+  void fromCoordAToCoordB(const Vec3r& p,
+			  Vec3r& q,
+			  const real transform[4][4]);
+
+  void fromWorldToCamera(const Vec3r& p,
+			 Vec3r& q,
+			 const real model_view_matrix[4][4]) {
+    fromCoordAToCoordB(p, q, model_view_matrix);
+  }
+
+  void fromCameraToRetina(const Vec3r& p,
+			  Vec3r& q,
+			  const real projection_matrix[4][4]) {
+    fromCoordAToCoordB(p, q, projection_matrix);
+  }
+
+  void fromRetinaToImage(const Vec3r& p,
+			 Vec3r& q,
+			 const int viewport[4]) {
+    // winX:
+    q[0] = viewport[0] + viewport[2] * (p[0] + 1.0) / 2.0;
+
+    // winY:
+    q[1] = viewport[1] + viewport[3] * (p[1] + 1.0) / 2.0;
+
+    // winZ:
+    q[2] = (p[2] + 1.0) / 2.0;
+  }
+
+  void fromWorldToImage(const Vec3r& p,
+			Vec3r& q,
+			const real model_view_matrix[4][4],
+			const real projection_matrix[4][4],
+			const int viewport[4]) {
+    Vec3r p1, p2;
+    fromWorldToCamera(p, p1, model_view_matrix);
+    fromCameraToRetina(p1, p2, projection_matrix);
+    fromRetinaToImage(p2, q, viewport);
+    q[2] = p1[2];
+  }
+
+  void fromWorldToImage(const Vec3r& p,
+			Vec3r& q,
+			const real transform[4][4],
+			const int viewport[4]) {
+    fromCoordAToCoordB(p, q, transform);
+
+    // winX:
+    q[0] = viewport[0] + viewport[2] * (q[0] + 1.0) / 2.0;
+
+    //winY:
+    q[1] = viewport[1] + viewport[3] * (q[1] + 1.0) / 2.0;
+  }
+
+  void fromImageToRetina(const Vec3r& p,
+			 Vec3r& q,
+			 const int viewport[4]) {
+    q = p;
+    q[0] = 2.0 * (q[0] - viewport[0]) / viewport[2] - 1;
+    q[1] = 2.0 * (q[1] - viewport[1]) / viewport[3] - 1;
+  }
+
+  void fromRetinaToCamera(const Vec3r& p,
+			  Vec3r& q,
+			  real z,
+			  const real projection_matrix[4][4]) {
+    q[0] = (-p[0] * z) / projection_matrix[0][0];
+    q[1] = (-p[1] * z) / projection_matrix[1][1];
+    q[2] = z;
+  }
+
+  void fromCameraToWorld(const Vec3r& p,
+			 Vec3r& q,
+			 const real model_view_matrix[4][4]) {
+
+    real translation[3] = { model_view_matrix[0][3],
+			      model_view_matrix[1][3],
+			      model_view_matrix[2][3] };
+    for (unsigned i = 0; i < 3; i++) {
+      q[i] = 0.0;
+      for (unsigned short j = 0; j < 3; j++)
+	q[i] += model_view_matrix[j][i] * (p[j] - translation[j]);
+    }
+  }
+
+
+  //
+  // Internal code
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  // Copyright 2001, softSurfer (www.softsurfer.com)
+  // This code may be freely used and modified for any purpose
+  // providing that this copyright notice is included with it.
+  // SoftSurfer makes no warranty for this code, and cannot be held
+  // liable for any real or imagined damage resulting from its use.
+  // Users of this code must verify correctness for their application.
+
+#define perp(u,v)  ((u)[0] * (v)[1] - (u)[1] * (v)[0])   // 2D perp product
+
+  inline bool intersect2dSegPoly(Vec2r* seg,
+				 Vec2r* poly,
+				 unsigned n) {
+    if (seg[0] == seg[1])
+      return false;
+
+    real  tE = 0;             // the maximum entering segment parameter
+    real  tL = 1;             // the minimum leaving segment parameter
+    real  t, N, D;            // intersect parameter t = N / D
+    Vec2r dseg;                // the segment direction vector
+    dseg = seg[1] - seg[0];
+    Vec2r e;                   // edge vector
+
+    for (unsigned i = 0; i < n; i++) { // process polygon edge poly[i]poly[i+1]
+      e = poly[i+1] - poly[i];
+      N = perp(e, seg[0] - poly[i]);
+      D = -perp(e, dseg);
+      if (fabs(D) < M_EPSILON) {
+	if (N < 0)
+	  return false;
+	else
+	  continue;
+      }
+
+      t = N / D;
+      if (D < 0) {       // segment seg is entering across this edge
+	if (t > tE) {    // new max tE
+	  tE = t;
+	  if (tE > tL)   // seg enters after leaving polygon
+	    return false;
+	}
+      }
+      else {             // segment seg is leaving across this edge
+	if (t < tL) {    // new min tL
+	  tL = t;
+	  if (tL < tE)   // seg leaves before entering polygon
+	    return false;
+	}
+      }
+    }
+
+    // tE <= tL implies that there is a valid intersection subsegment
+    return true;
+  }
+
+  inline bool overlapPlaneBox(Vec3r& normal, real d, Vec3r& maxbox) {
+    Vec3r vmin, vmax;
+
+    for(unsigned q = X; q <= Z; q++) {
+	if(normal[q] > 0.0f) {
+	  vmin[q] = -maxbox[q];
+	  vmax[q] = maxbox[q];
+	}
+	else {
+	  vmin[q] = maxbox[q];
+	  vmax[q] = -maxbox[q];
+	}
+    }
+    if((normal * vmin) + d > 0.0f)
+      return false;
+    if((normal * vmax) + d >= 0.0f)
+      return true;
+    return false;
+  }
+
+  inline void fromCoordAToCoordB(const Vec3r&p,
+				 Vec3r& q,
+				 const real transform[4][4]) {
+    HVec3r hp(p);
+    HVec3r hq(0, 0, 0, 0);
+
+    for (unsigned i = 0; i < 4; i++)
+      for (unsigned j = 0; j < 4; j++)
+	hq[i] += transform[i][j] * hp[j];
+
+    if(hq[3] == 0) {
+      q = p;
+      return;
+    }
+
+    for (unsigned k = 0; k < 3; k++)
+      q[k] = hq[k] / hq[3];
+  }
+
+} // end of namespace GeomUtils
diff --git a/source/blender/freestyle/intern/geometry/GeomUtils.h b/source/blender/freestyle/intern/geometry/GeomUtils.h
new file mode 100755
index 00000000000..53c94c22f8b
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/GeomUtils.h
@@ -0,0 +1,309 @@
+//
+//  Filename         : GeomUtils.h
+//  Author(s)        : Stephane Grabli
+//  Purpose          : Various tools for geometry
+//  Date of creation : 12/04/2002
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef  GEOMUTILS_H
+# define GEOMUTILS_H
+
+# include <vector>
+# include "../system/FreestyleConfig.h"
+# include "Geom.h"
+
+using namespace std;
+using namespace Geometry;
+
+namespace GeomUtils {
+
+  //
+  // Templated procedures
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  /*! Computes the distance from a point P to a segment AB */
+  template<class T>
+  real distPointSegment( const T& P, const T& A , const T& B) { 
+    T AB, AP, BP;
+    AB = B - A;
+    AP = P - A;
+    BP = P - B;
+
+    real c1(AB * AP);
+    if (c1 <= 0)
+      return AP.norm();
+
+    real c2(AB * AB);
+    if (c2 <= c1)
+      return BP.norm();
+
+    real b = c1 / c2;
+    T Pb, PPb;
+    Pb = A + b * AB;
+    PPb = P - Pb;
+
+    return PPb.norm();
+  } 
+
+  //
+  // Non-templated procedures
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  typedef enum {
+    DONT_INTERSECT,
+    DO_INTERSECT,
+    COLINEAR,
+    COINCIDENT
+  } intersection_test;
+
+  LIB_GEOMETRY_EXPORT
+  intersection_test intersect2dSeg2dSeg(const Vec2r& p1, const Vec2r& p2, // first segment
+					const Vec2r& p3, const Vec2r& p4, // second segment
+					Vec2r& res);                      // found intersection point
+
+  LIB_GEOMETRY_EXPORT
+  intersection_test intersect2dLine2dLine(const Vec2r& p1, const Vec2r& p2, // first segment
+					  const Vec2r& p3, const Vec2r& p4, // second segment
+					  Vec2r& res);                      // found intersection point
+
+  LIB_GEOMETRY_EXPORT
+  intersection_test intersect2dSeg2dSegParametric(const Vec2r& p1, const Vec2r& p2, // first segment
+						  const Vec2r& p3, const Vec2r& p4, // second segment
+						  real& t,                          // I = P1 + t * P1P2)
+						  real& u);                         // I = P3 + u * P3P4
+
+  /*! check whether a 2D segment intersect a 2D region or not */
+  LIB_GEOMETRY_EXPORT
+  bool intersect2dSeg2dArea(const Vec2r& min,
+			    const Vec2r& max,
+			    const Vec2r& A,
+			    const Vec2r& B);
+  
+  /*! check whether a 2D segment is included in a 2D region or not */
+  LIB_GEOMETRY_EXPORT
+  bool include2dSeg2dArea(const Vec2r& min,
+			  const Vec2r& max,
+			  const Vec2r& A,
+			  const Vec2r& B);
+
+  /*! Box-triangle overlap test, adapted from Tomas Akenine-Möller code */
+  LIB_GEOMETRY_EXPORT
+  bool overlapTriangleBox(Vec3r& boxcenter,
+			  Vec3r& boxhalfsize,
+			  Vec3r triverts[3]);
+
+  /*! Fast, Minimum Storage Ray-Triangle Intersection,
+   * adapted from Tomas Möller and Ben Trumbore code.
+   */
+  LIB_GEOMETRY_EXPORT
+  bool intersectRayTriangle(Vec3r& orig, Vec3r& dir,
+			    Vec3r& v0, Vec3r& v1, Vec3r& v2,
+			    real& t,                   // I = orig + t * dir
+                real& u, real& v,          // I = (1-u-v)*v0+u*v1+v*v2
+			    real epsilon = M_EPSILON); // the epsilon to use
+
+  /*! Intersection between plane and ray
+   * adapted from Graphics Gems, Didier Badouel 
+   */
+  LIB_GEOMETRY_EXPORT
+  intersection_test intersectRayPlane(Vec3r& orig, Vec3r& dir,   // ray origin and direction
+				      Vec3r& norm, real d,       // plane's normal and offset (plane = { P / P.N + d = 0 })
+				      real& t,                   // I = orig + t * dir
+				      real epsilon = M_EPSILON); // the epsilon to use
+
+ /*! Intersection Ray-Bounding box (axis aligned).
+  *  Adapted from Williams et al, "An Efficient Robust Ray-Box Intersection Algorithm", 
+  *  JGT 10:1 (2005), pp. 49-54.
+  *  Returns
+  */
+  LIB_GEOMETRY_EXPORT
+      bool intersectRayBBox(const Vec3r& orig, const Vec3r& dir,   // ray origin and direction
+      const Vec3r& boxMin, const Vec3r& boxMax, // the bbox
+      real t0, real t1,     // the interval in which at least on of the intersections must happen
+      real& tmin, real& tmax,                  // Imin=orig+tmin*dir is the first intersection, Imax=orig+tmax*dir is the second intersection
+      real epsilon = M_EPSILON); // the epsilon to use
+
+
+  /*! Checks whether 3D point P lies inside or outside of the triangle ABC */
+  LIB_GEOMETRY_EXPORT
+  bool includePointTriangle(Vec3r& P,
+			    Vec3r& A,
+			    Vec3r& B,
+			    Vec3r& C);
+
+  LIB_GEOMETRY_EXPORT
+  void transformVertex(const Vec3r& vert,
+		       const Matrix44r& matrix,
+		       Vec3r& res);
+  
+  LIB_GEOMETRY_EXPORT
+  void transformVertices(const vector<Vec3r>& vertices,
+			 const Matrix44r& trans,
+			 vector<Vec3r>& res);
+
+  LIB_GEOMETRY_EXPORT
+  Vec3r rotateVector(const Matrix44r& mat, const Vec3r& v);
+
+  //
+  // Coordinates systems changing procedures
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  /*! From world to image
+   *  p
+   *    point's coordinates expressed in world coordinates system
+   *  q
+   *    vector in which the result will be stored
+   *  model_view_matrix
+   *    The model view matrix expressed in line major order (OpenGL
+   *    matrices are column major ordered)
+   *  projection_matrix
+   *    The projection matrix expressed in line major order (OpenGL
+   *    matrices are column major ordered)
+   *  viewport
+   *    The viewport: x,y coordinates followed by width and height (OpenGL like viewport)
+   */
+  LIB_GEOMETRY_EXPORT
+  void fromWorldToImage(const Vec3r& p,
+			Vec3r& q,
+			const real model_view_matrix[4][4],
+			const real projection_matrix[4][4],
+			const int viewport[4]);
+
+  /*! From world to image
+   *  p
+   *    point's coordinates expressed in world coordinates system
+   *  q
+   *    vector in which the result will be stored
+   *  transform
+   *    The transformation matrix (gathering model view and projection), 
+   *    expressed in line major order (OpenGL matrices are column major ordered)
+   *  viewport
+   *    The viewport: x,y coordinates followed by width and height (OpenGL like viewport)
+   */
+  LIB_GEOMETRY_EXPORT
+  void fromWorldToImage(const Vec3r& p,
+			Vec3r& q,
+			const real transform[4][4],
+			const int viewport[4]);
+
+  /*! Projects from world coordinates to camera coordinates 
+   *  Returns the point's coordinates expressed in the camera's 
+   *  coordinates system.
+   *  p
+   *    point's coordinates expressed in world coordinates system
+   *  q
+   *    vector in which the result will be stored
+   *  model_view_matrix
+   *    The model view matrix expressed in line major order (OpenGL
+   *    matrices are column major ordered)
+   */
+  LIB_GEOMETRY_EXPORT
+  void fromWorldToCamera(const Vec3r& p,
+			 Vec3r& q,
+			 const real model_view_matrix[4][4]);
+  
+  /*! Projects from World Coordinates to retina coordinates
+   *  Returns the point's coordinates expressed in Retina system.
+   *  p
+   *    point's coordinates expressed in camera system
+   *  q
+   *    vector in which the result will be stored
+   *  projection_matrix
+   *    The projection matrix expressed in line major order (OpenGL
+   *    matrices are column major ordered)
+   */
+  LIB_GEOMETRY_EXPORT
+  void fromCameraToRetina(const Vec3r& p,
+			  Vec3r& q,
+			  const real projection_matrix[4][4]);
+  
+  /*! From retina to image.
+   *  Returns the coordinates expressed in Image coorinates system.
+   *  p
+   *    point's coordinates expressed in retina system
+   *  q
+   *    vector in which the result will be stored
+   *  viewport
+   *    The viewport: x,y coordinates followed by width and height (OpenGL like viewport).
+   */
+  LIB_GEOMETRY_EXPORT
+  void fromRetinaToImage(const Vec3r& p,
+			 Vec3r& q,
+			 const int viewport[4]);
+  
+  /*! From image to retina
+   *  p
+   *    point's coordinates expressed in image system
+   *  q
+   *    vector in which the result will be stored
+   *  viewport
+   *    The viewport: x,y coordinates followed by width and height (OpenGL like viewport).
+   */
+  LIB_GEOMETRY_EXPORT
+  void fromImageToRetina(const Vec3r& p,
+			 Vec3r& q,
+			 const int viewport[4]);
+ 
+  /*! computes the coordinates of q in the camera coordinates system, 
+   *  using the known z coordinates of the 3D point. 
+   *  That means that this method does not inverse any matrices, 
+   *  it only computes X and Y from x,y and Z)
+   *  p
+   *    point's coordinates expressed in retina system
+   *  q
+   *    vector in which the result will be stored
+   *  projection_matrix
+   *    The projection matrix expressed in line major order (OpenGL
+   *    matrices are column major ordered)
+   
+   */
+  LIB_GEOMETRY_EXPORT
+  void fromRetinaToCamera(const Vec3r& p,
+			  Vec3r& q,
+			  real z,
+			  const real projection_matrix[4][4]);
+
+  /*! Projects from camera coordinates to world coordinates
+   *  Returns the point's coordinates expressed in the world's 
+   *  coordinates system.
+   *  p
+   *    point's coordinates expressed in the camera coordinates system
+   *  q
+   *    vector in which the result will be stored
+   *  model_view_matrix
+   *    The model view matrix expressed in line major order (OpenGL
+   *    matrices are column major ordered)
+   */
+  LIB_GEOMETRY_EXPORT
+  void fromCameraToWorld(const Vec3r& p,
+			 Vec3r& q,
+			 const real model_view_matrix[4][4]);
+
+} // end of namespace GeomUtils
+
+#endif // GEOMUTILS_H
diff --git a/source/blender/freestyle/intern/geometry/Grid.cpp b/source/blender/freestyle/intern/geometry/Grid.cpp
new file mode 100755
index 00000000000..59b730358bc
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/Grid.cpp
@@ -0,0 +1,388 @@
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#include "Grid.h"
+#include "BBox.h"
+#include <cassert>
+#include <stdexcept>
+
+// Grid Visitors
+/////////////////
+void allOccludersGridVisitor::examineOccluder(Polygon3r *occ){
+    occluders_.push_back(occ);
+}
+
+bool inBox(const Vec3r& inter, const Vec3r& box_min, const Vec3r& box_max){
+    if(((inter.x()>=box_min.x()) && (inter.x() <box_max.x()))
+        && ((inter.y()>=box_min.y()) && (inter.y() <box_max.y()))
+        && ((inter.z()>=box_min.z()) && (inter.z() <box_max.z()))
+        ){
+            return true;
+        }
+        return false;
+}
+void firstIntersectionGridVisitor::examineOccluder(Polygon3r *occ){
+
+    // check whether the edge and the polygon plane are coincident:
+    //-------------------------------------------------------------
+    //first let us compute the plane equation.
+    Vec3r v1(((occ)->getVertices())[0]);           
+    Vec3d normal((occ)->getNormal());
+    double d = -(v1 * normal);
+
+    double tmp_u, tmp_v, tmp_t;
+    if((occ)->rayIntersect(ray_org_, ray_dir_, tmp_t, tmp_u, tmp_v)){
+        if (fabs(ray_dir_ * normal) > 0.0001){
+            // Check whether the intersection is in the cell:
+            if(inBox(ray_org_+tmp_t*ray_dir_/ray_dir_.norm(), current_cell_->getOrigin(), current_cell_->getOrigin()+cell_size_)){
+
+                //Vec3d bboxdiag(_scene3d->bbox().getMax()-_scene3d->bbox().getMin());
+                //if ((t>1.0E-06*(min(min(bboxdiag.x(),bboxdiag.y()),bboxdiag.z()))) && (t<raylength)){
+                if(tmp_t < t_){
+                    occluder_ = occ;
+                    u_ = tmp_u;
+                    v_ = tmp_v;
+                    t_ = tmp_t;
+                }
+            }else{
+                occ->userdata2 = 0;
+            }
+        }
+    }
+}
+
+bool firstIntersectionGridVisitor::stop(){
+    if(occluder_)
+        return true;
+    return false;
+}
+
+// Grid
+/////////////////
+
+void Grid::clear() {
+  if (_occluders.size() != 0) {
+    for(OccludersSet::iterator it = _occluders.begin();
+	it != _occluders.end();
+	it++) {
+      delete (*it);
+    }
+    _occluders.clear();
+  }
+
+  _size = Vec3r(0, 0, 0);
+  _cell_size = Vec3r(0, 0, 0);
+  _orig = Vec3r(0, 0, 0);
+  _cells_nb = Vec3u(0, 0, 0);
+  //_ray_occluders.clear();
+}       
+
+void Grid::configure(const Vec3r& orig,
+		     const Vec3r& size,
+		     unsigned nb) {
+                 
+  _orig = orig;
+ Vec3r tmpSize=size;
+  // Compute the volume of the desired grid
+  real grid_vol = size[0] * size[1] * size[2];
+
+  if(grid_vol == 0){
+    double min=DBL_MAX;
+    int index;
+    int nzeros=0;
+    for(int i=0;i<3;++i){
+        if(size[i] == 0){
+            ++nzeros;
+            index=i;
+        }
+        if((size[i]!=0) && (min>size[i])){
+            min=size[i];
+        }
+    }
+    if(nzeros>1){
+        throw std::runtime_error("Warning: the 3D grid has more than one null dimension");
+    }
+    tmpSize[index]=min;
+    _orig[index] = _orig[index]-min/2;
+  }
+  // Compute the desired volume of a single cell
+  real cell_vol = grid_vol / nb;
+  // The edge of such a cubic cell is cubic root of cellVolume
+  real edge = pow(cell_vol, 1.0 / 3.0);
+
+  // We compute the number of cells par edge
+  // such as we cover at least the whole box.
+  unsigned i;
+  for (i = 0; i < 3; i++)
+    _cells_nb[i] = (unsigned)floor(tmpSize[i] / edge) + 1;
+
+  _size = tmpSize;
+
+  for(i = 0; i < 3; i++)
+    _cell_size[i] = _size[i] / _cells_nb[i];
+}
+
+void Grid::insertOccluder(Polygon3r* occluder) {
+  const vector<Vec3r> vertices = occluder->getVertices();
+  if (vertices.size() == 0)
+    return;
+
+  // add this occluder to the grid's occluders list
+  addOccluder(occluder);
+
+  // find the bbox associated to this polygon
+  Vec3r min, max;
+  occluder->getBBox(min, max);
+
+  // Retrieve the cell x, y, z cordinates associated with these min and max
+  Vec3u imax, imin;
+  getCellCoordinates(max, imax);
+  getCellCoordinates(min, imin);
+  
+  // We are now going to fill in the cells overlapping with the
+  // polygon bbox.
+  // If the polygon is a triangle (most of cases), we also
+  // check for each of these cells if it is overlapping with
+  // the triangle in order to only fill in the ones really overlapping
+  // the triangle.
+
+  unsigned i, x, y, z;
+  vector<Vec3r>::const_iterator it;
+  Vec3u coord;
+
+  if (vertices.size() == 3) { // Triangle case
+    Vec3r triverts[3];
+    i = 0;
+    for(it = vertices.begin();
+	it != vertices.end();
+	it++) {
+      triverts[i] = Vec3r(*it);
+      i++;
+    }
+
+    Vec3r boxmin, boxmax;
+
+    for (z = imin[2]; z <= imax[2]; z++)
+      for (y = imin[1]; y <= imax[1]; y++)
+	for (x = imin[0]; x <= imax[0]; x++) {
+	  coord[0] = x;
+	  coord[1] = y;
+	  coord[2] = z;
+	  // We retrieve the box coordinates of the current cell
+	  getCellBox(coord, boxmin, boxmax);
+	  // We check whether the triangle and the box ovewrlap:
+	  Vec3r boxcenter((boxmin + boxmax) / 2.0);
+	  Vec3r boxhalfsize(_cell_size / 2.0);
+	  if (GeomUtils::overlapTriangleBox(boxcenter, boxhalfsize, triverts)) {
+	    // We must then create the Cell and add it to the cells list
+	    // if it does not exist yet.
+	    // We must then add the occluder to the occluders list of this cell.
+	    Cell* cell = getCell(coord);
+	    if (!cell) {
+	      cell = new Cell(boxmin);
+	      fillCell(coord, *cell);
+	    }
+	    cell->addOccluder(occluder);
+	  }
+	}
+  }
+  else { // The polygon is not a triangle, we add all the cells overlapping the polygon bbox.
+    for (z = imin[2]; z <= imax[2]; z++)
+      for (y = imin[1]; y <= imax[1]; y++)
+	for (x = imin[0]; x <= imax[0]; x++) {
+	  coord[0] = x;
+	  coord[1] = y;
+	  coord[2] = z;
+	  Cell* cell = getCell(coord);
+	  if (!cell) {
+	    Vec3r orig;
+	    getCellOrigin(coord, orig);
+	    cell = new Cell(orig);
+	    fillCell(coord, *cell);
+	  }
+	  cell->addOccluder(occluder);
+	}
+  }
+}
+
+bool Grid::nextRayCell(Vec3u& current_cell, Vec3u& next_cell) {
+  next_cell = current_cell;
+  real t_min, t;
+  unsigned i;
+ 
+  t_min = FLT_MAX;    // init tmin with handle of the case where one or 2 _u[i] = 0. 
+  unsigned coord = 0; // predominant coord(0=x, 1=y, 2=z)
+
+
+  // using a parametric equation of
+  // a line : B = A + t u, we find
+  // the tx, ty and tz respectively coresponding
+  // to the intersections with the plans:
+  // x = _cell_size[0], y = _cell_size[1], z = _cell_size[2]
+  for (i = 0; i < 3; i++) {
+    if (_ray_dir[i] == 0)
+      continue;
+    if (_ray_dir[i] > 0)
+      t = (_cell_size[i] - _pt[i]) / _ray_dir[i];
+    else
+      t = -_pt[i] / _ray_dir[i];
+    if (t < t_min) {
+      t_min = t;
+      coord = i;
+    }
+  }
+
+  // We use the parametric line equation and
+  // the found t (tamx) to compute the
+  // B coordinates:
+  Vec3r pt_tmp(_pt);
+  _pt = pt_tmp + t_min * _ray_dir;
+    
+  // We express B coordinates in the next cell
+  // coordinates system. We just have to
+  // set the coordinate coord of B to 0
+  // of _CellSize[coord] depending on the sign
+  // of _u[coord]
+  if (_ray_dir[coord] > 0) {
+    next_cell[coord]++;
+    _pt[coord] -= _cell_size[coord];
+    // if we are out of the grid, we must stop
+    if (next_cell[coord] >= _cells_nb[coord])
+      return false;
+  }
+  else {
+    int tmp = next_cell[coord] - 1;
+    _pt[coord] = _cell_size[coord];
+    if (tmp < 0)
+      return false;
+    next_cell[coord]--;
+  }
+
+  _t += t_min;
+  if (_t >= _t_end)
+    return false;
+
+  return true;
+}
+
+void Grid::castRay(const Vec3r& orig,
+		   const Vec3r& end,
+		   OccludersSet& occluders,
+		   unsigned timestamp) {
+  initRay(orig, end, timestamp);
+  allOccludersGridVisitor visitor(occluders);
+  castRayInternal(visitor);
+}
+
+void Grid::castInfiniteRay(const Vec3r& orig,
+			   const Vec3r& dir,
+			   OccludersSet& occluders,
+			   unsigned timestamp) {
+  Vec3r end = Vec3r(orig + FLT_MAX * dir / dir.norm());
+  bool inter = initInfiniteRay(orig, dir, timestamp);
+  if(!inter)
+      return;
+  allOccludersGridVisitor visitor(occluders);
+  castRayInternal(visitor);
+}
+  
+Polygon3r* Grid::castRayToFindFirstIntersection(const Vec3r& orig,
+                   const Vec3r& dir,
+                   double& t,
+                   double& u,
+                   double& v,
+                   unsigned timestamp){
+    Polygon3r *occluder = 0;
+    Vec3r end = Vec3r(orig + FLT_MAX * dir / dir.norm());
+    bool inter = initInfiniteRay(orig, dir, timestamp);
+    if(!inter){
+        return 0;
+    }
+    firstIntersectionGridVisitor visitor(orig,dir,_cell_size);
+    castRayInternal(visitor);
+    occluder = visitor.occluder();
+    t = visitor.t_;
+    u = visitor.u_;
+    v = visitor.v_;
+    return occluder;
+}
+
+void Grid::initRay (const Vec3r &orig,
+		    const Vec3r& end,
+		    unsigned timestamp) {
+  _ray_dir = end - orig;
+  _t_end = _ray_dir.norm();
+  _t = 0;
+  _ray_dir.normalize();
+  _timestamp = timestamp;
+
+  for(unsigned i = 0; i < 3; i++) {
+    _current_cell[i] = (unsigned)floor((orig[i] - _orig[i]) / _cell_size[i]);
+    unsigned u = _current_cell[i];
+    _pt[i] = orig[i] - _orig[i] - _current_cell[i] * _cell_size[i];
+  }
+  //_ray_occluders.clear();
+
+}
+
+bool Grid::initInfiniteRay (const Vec3r &orig,
+		    const Vec3r& dir,
+		    unsigned timestamp) {
+  _ray_dir = dir;
+  _t_end = FLT_MAX;
+  _t = 0;
+  _ray_dir.normalize();
+  _timestamp = timestamp;
+
+  // check whether the origin is in or out the box:
+  Vec3r boxMin(_orig);
+  Vec3r boxMax(_orig+_size);
+  BBox<Vec3r> box(boxMin, boxMax);
+  if(box.inside(orig)){
+      for(unsigned i = 0; i < 3; i++) {
+          _current_cell[i] = (unsigned)floor((orig[i] - _orig[i]) / _cell_size[i]);
+          unsigned u = _current_cell[i];
+          _pt[i] = orig[i] - _orig[i] - _current_cell[i] * _cell_size[i];
+      }
+  }else{
+      // is the ray intersecting the box?
+      real tmin(-1.0), tmax(-1.0);
+      if(GeomUtils::intersectRayBBox(orig, _ray_dir, boxMin, boxMax, 0, _t_end, tmin, tmax)){
+        assert(tmin != -1.0);
+        Vec3r newOrig = orig + tmin*_ray_dir;
+        for(unsigned i = 0; i < 3; i++) {
+            _current_cell[i] = (unsigned)floor((newOrig[i] - _orig[i]) / _cell_size[i]);
+            if(_current_cell[i] == _cells_nb[i])
+                _current_cell[i] = _cells_nb[i] - 1;
+            unsigned u = _current_cell[i];
+            _pt[i] = newOrig[i] - _orig[i] - _current_cell[i] * _cell_size[i];
+        }
+
+      }else{
+          return false;
+      }
+  }
+  //_ray_occluders.clear();
+
+  return true;
+
+}
+
diff --git a/source/blender/freestyle/intern/geometry/Grid.h b/source/blender/freestyle/intern/geometry/Grid.h
new file mode 100755
index 00000000000..6197721bb45
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/Grid.h
@@ -0,0 +1,358 @@
+//
+//  Filename         : Grid.h
+//  Author(s)        : Stephane Grabli
+//  Purpose          : Base class to define a cell grid surrounding
+//                     the bounding box of the scene
+//  Date of creation : 30/07/2002
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef  GRID_H
+# define GRID_H
+
+# include <float.h>
+# include <vector>
+# include "../system/FreestyleConfig.h"
+# include "GeomUtils.h"
+# include "Geom.h"
+# include "Polygon.h"
+
+using namespace std;
+using namespace Geometry;
+
+typedef vector<Polygon3r*>	OccludersSet;
+
+
+//
+// Class to define cells used by the regular grid
+//
+///////////////////////////////////////////////////////////////////////////////
+
+class LIB_GEOMETRY_EXPORT Cell
+{
+ public:
+  
+  Cell(Vec3r& orig) {
+    _orig = orig;
+  }
+
+  virtual ~Cell() {}
+
+  inline void addOccluder(Polygon3r* o) {
+    if (o)
+      _occluders.push_back(o);
+  }
+
+  inline const Vec3r& getOrigin() {
+    return _orig;
+  }
+
+  inline OccludersSet& getOccluders() {
+    return _occluders;
+  }
+  
+ private:
+
+  Vec3r		_orig;
+  OccludersSet	_occluders;
+};
+
+
+class GridVisitor{
+public:
+    virtual void discoverCell(Cell *cell) {}
+    virtual void examineOccluder(Polygon3r *occ) {}
+    virtual void finishCell(Cell *cell) {}
+    virtual bool stop() {return false;}
+};
+
+/*! Gathers all the occluders belonging to the cells
+ * traversed by the ray */
+class allOccludersGridVisitor : public GridVisitor{
+public:
+    allOccludersGridVisitor(OccludersSet& occluders)
+        :GridVisitor(), occluders_(occluders){}
+    virtual void examineOccluder(Polygon3r *occ);
+
+    OccludersSet& occluders() {return occluders_;}
+    void clear() {occluders_.clear();}
+
+private:
+    OccludersSet& occluders_;
+};
+
+/*! Finds the first intersection and breaks. The occluder and 
+ * the intersection information are stored and accessible.
+ */
+class firstIntersectionGridVisitor : public GridVisitor {
+public:
+      firstIntersectionGridVisitor(const Vec3r& ray_org, const Vec3r& ray_dir, const Vec3r& cell_size) : 
+      GridVisitor(), ray_org_(ray_org), cell_size_(cell_size),ray_dir_(ray_dir),occluder_(0),
+            u_(0),v_(0),t_(DBL_MAX),current_cell_(0){}
+      virtual ~firstIntersectionGridVisitor() {}
+
+    virtual void discoverCell(Cell *cell) {current_cell_=cell;}
+    virtual void examineOccluder(Polygon3r *occ);
+    virtual bool stop();
+
+    Polygon3r * occluder() {return occluder_;}
+    
+public:
+    double u_, v_, t_;
+private:
+    Polygon3r *occluder_;
+    Vec3r ray_org_, ray_dir_;
+    Vec3r cell_size_;
+    Cell * current_cell_;
+};
+
+//
+// Class to define a regular grid used for ray casting computations
+//
+///////////////////////////////////////////////////////////////////////////////
+
+class LIB_GEOMETRY_EXPORT Grid
+{
+public:
+
+  /*! Builds a Grid 
+   *  Must be followed by a call to configure()
+   */
+  Grid() {}
+
+  virtual ~Grid() {
+    clear();
+  }
+
+  /*! clears the grid
+   *  Deletes all the cells, clears the hashtable,
+   *  resets size, size of cell, number of cells.
+   */
+  virtual void clear();
+
+  /*! Sets the different parameters of the grid
+   *    orig
+   *      The grid origin
+   *    size
+   *      The grid's dimensions
+   *    nb
+   *      The number of cells of the grid
+   */
+  virtual void configure(const Vec3r& orig, const Vec3r& size, unsigned nb);
+  
+  /*! returns a vector of integer containing the
+   *  coordinates of the cell containing the point
+   *  passed as argument
+   *    p
+   *      The point for which we're looking the cell
+   */
+  inline void getCellCoordinates(const Vec3r& p, Vec3u& res) {
+    int tmp;
+    for (int i = 0; i < 3; i++) {
+      tmp = (int)((p[i] - _orig[i]) / _cell_size[i]);
+      if (tmp < 0)
+	res[i] = 0;
+      else if ((unsigned)tmp >= _cells_nb[i])
+	res[i] = _cells_nb[i] - 1;
+      else
+	res[i] = tmp;
+    }
+  }
+
+  /*! Fills the case corresponding to coord with the cell */
+  virtual void fillCell(const Vec3u& coord, Cell& cell) = 0;
+
+  /*! returns the cell whose coordinates
+   *  are pased as argument
+   */
+  virtual Cell* getCell(const Vec3u& coord) = 0;
+
+  /*! returns the cell containing the point
+   *  passed as argument. If the cell is empty
+   *  (contains no occluder), NULL is returned
+   *    p
+   *      The point for which we're looking the cell
+   */
+  inline Cell* getCell(const Vec3r& p) {    
+    Vec3u coord;
+    getCellCoordinates(p, coord);
+    return getCell(coord);
+  }
+
+  /*! Retrieves the x,y,z coordinates of the origin of the cell whose coordinates (i,j,k)
+   *  is passed as argument
+   *    cell_coord
+   *      i,j,k integer coordinates for the cell
+   *    orig
+   *      x,y,x vector to be filled in with the cell origin's coordinates
+   */
+  inline void getCellOrigin(const Vec3u& cell_coord, Vec3r& orig) {
+    for (unsigned i = 0; i < 3; i++)
+      orig[i] = _orig[i] + cell_coord[i] * _cell_size[i];
+  }
+
+  /*! Retrieves the box corresponding to the cell whose coordinates
+   *  are passed as argument.
+   *    cell_coord
+   *      i,j,k integer coordinates for the cell
+   *    min_out
+   *      The min x,y,x vector of the box. Filled in by the method.
+   *    max_out
+   *      The max x,y,z coordinates of the box. Filled in by the method.
+   */
+  inline void getCellBox(const Vec3u& cell_coord, Vec3r& min_out, Vec3r& max_out) {
+    getCellOrigin(cell_coord, min_out);
+    max_out = min_out + _cell_size;
+  }
+
+  /*! inserts a convex polygon occluder
+   *  This method is quite coarse insofar as it
+   *  adds all cells intersecting the polygon bounding box
+   *    convex_poly
+   *      The list of 3D points constituing a convex polygon
+   */
+  void insertOccluder(Polygon3r * convex_poly);
+
+  /*! Adds an occluder to the list of occluders */
+  void addOccluder(Polygon3r* occluder) {
+    _occluders.push_back(occluder);
+  }
+
+  /*! Casts a ray between a starting point and an ending point
+   *  Returns the list of occluders contained
+   *  in the cells intersected by this ray
+   *  Starts with a call to InitRay.
+   */
+  void castRay(const Vec3r& orig,
+	       const Vec3r& end,
+	       OccludersSet& occluders,
+	       unsigned timestamp);
+
+  /*! Casts an infinite ray (still finishing at the end of the grid) from a starting point and in a given direction.
+   *  Returns the list of occluders contained
+   *  in the cells intersected by this ray
+   *  Starts with a call to InitRay.
+   */
+  void castInfiniteRay(const Vec3r& orig,
+		       const Vec3r& dir,
+		       OccludersSet& occluders,
+		       unsigned timestamp);
+
+  /*! Casts an infinite ray (still finishing at the end of the grid) from a starting point and in a given direction.
+  *  Returns the first intersection (occluder,t,u,v) or null.
+  *  Starts with a call to InitRay.
+  */
+  Polygon3r * castRayToFindFirstIntersection(const Vec3r& orig,
+      const Vec3r& dir,
+      double& t,
+      double& u,
+      double& v,
+      unsigned timestamp);
+
+
+  /*! Init all structures and values for computing
+   *  the cells intersected by this new ray
+   */
+  void initRay (const Vec3r &orig,
+		const Vec3r& end,
+		unsigned timestamp);
+
+  /*! Init all structures and values for computing
+   *  the cells intersected by this infinite ray.
+   * Returns false if the ray doesn't intersect the 
+   * grid.
+   */
+  bool initInfiniteRay (const Vec3r &orig,
+		const Vec3r& dir,
+		unsigned timestamp);
+
+ 
+  /*! Accessors */
+  inline const Vec3r& getOrigin() const {
+    return _orig;
+  }
+  inline Vec3r gridSize() const {
+    return _size;
+  }
+  inline Vec3r getCellSize() const {
+    return _cell_size;
+  }
+
+  void displayDebug() {
+    cerr << "Cells nb     : " << _cells_nb << endl;
+    cerr << "Cell size    : " << _cell_size << endl;
+    cerr << "Origin       : " << _orig << endl;
+    cerr << "Occluders nb : " << _occluders.size() << endl;
+  }
+
+ protected:
+
+  /*! Core of castRay and castInfiniteRay, find occluders
+   *  along the given ray
+   */
+  inline void castRayInternal(GridVisitor& visitor) {
+    Cell* current_cell = NULL;
+    do {
+      current_cell = getCell(_current_cell);
+      if (current_cell){
+          visitor.discoverCell(current_cell);
+          OccludersSet& occluders = current_cell->getOccluders(); // FIXME: I had forgotten the ref &
+          for (OccludersSet::iterator it = occluders.begin();
+              it != occluders.end();
+              it++) {
+                  if ((unsigned)(*it)->userdata2 != _timestamp) {
+                      (*it)->userdata2 = (void*)_timestamp;
+                      visitor.examineOccluder(*it);
+                  }
+              }
+          visitor.finishCell(current_cell);
+      }
+    } while ((!visitor.stop()) && (nextRayCell(_current_cell, _current_cell)));
+  }
+
+ 
+  /*! returns the  cell next to the cell
+   *  passed as argument.
+   */
+  bool nextRayCell(Vec3u& current_cell, Vec3u& next_cell);
+
+  unsigned	_timestamp;
+
+  Vec3u		_cells_nb;  // number of cells for x,y,z axis
+  Vec3r		_cell_size; // cell x,y,z dimensions
+  Vec3r		_size;      // grid x,y,x dimensions
+  Vec3r		_orig;      // grid origin
+
+  Vec3r		_ray_dir;      // direction vector for the ray
+  Vec3u		_current_cell; // The current cell being processed (designated by its 3 coordinates)
+  Vec3r		_pt;           // Points corresponding to the incoming and outgoing intersections
+                               // of one cell with the ray
+  real		_t_end;        // To know when we are at the end of the ray
+  real		_t;
+
+  //OccludersSet _ray_occluders; // Set storing the occluders contained in the cells traversed by a ray
+  OccludersSet _occluders;     // List of all occluders inserted in the grid
+};
+
+#endif // GRID_H
diff --git a/source/blender/freestyle/intern/geometry/HashGrid.cpp b/source/blender/freestyle/intern/geometry/HashGrid.cpp
new file mode 100755
index 00000000000..3cf845d57ef
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/HashGrid.cpp
@@ -0,0 +1,41 @@
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#include "HashGrid.h"
+
+void HashGrid::clear()
+{
+  if(!_cells.empty()) {
+    for(GridHashTable::iterator it = _cells.begin();
+	it !=_cells.end();
+	it++) {
+      Cell* cell = (*it).second;
+      delete cell;
+    }
+    _cells.clear();
+  }
+
+  Grid::clear();
+}
+
+void HashGrid::configure(const Vec3r& orig, const Vec3r& size, unsigned nb) {
+  Grid::configure(orig, size, nb);
+}
diff --git a/source/blender/freestyle/intern/geometry/HashGrid.h b/source/blender/freestyle/intern/geometry/HashGrid.h
new file mode 100755
index 00000000000..f6605957676
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/HashGrid.h
@@ -0,0 +1,109 @@
+//
+//  Filename         : HashGrid.h
+//  Author(s)        : Stephane Grabli
+//  Purpose          : Class to define a cell grid surrounding the
+//                     bounding box of the scene
+//  Date of creation : 30/07/2002
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef  HASHGRID_H
+# define HASHGRID_H
+
+//# if defined(__GNUC__) && (__GNUC__ >= 3)
+// hash_map is not part of the C++ standard anymore;
+// hash_map.h has been kept though for backward compatibility
+//#  include <hash_map.h>
+//# else
+//#  include <hash_map>
+//# endif
+
+# include "Grid.h"
+# include <map>
+/*! Defines a hash table used for searching the Cells */
+struct GridHasher{
+#define _MUL 950706376UL
+#define _MOD 2147483647UL
+  inline size_t operator() (const Vec3u& p) const 
+  {
+    size_t res = ((unsigned long) (p[0] * _MUL)) % _MOD;
+    res = ((res + (unsigned long) (p[1]) * _MUL)) % _MOD;
+    return ((res +(unsigned long) (p[2]) * _MUL)) % _MOD;
+  }
+};
+
+/*! Class to define a regular grid used for ray
+  casting computations */
+
+class LIB_GEOMETRY_EXPORT HashGrid : public Grid
+{
+ public:
+
+  typedef map<Vec3u, Cell*> GridHashTable;
+  
+  HashGrid() : Grid() {}
+
+  virtual ~HashGrid() {
+    clear();
+  }
+
+  /*! clears the grid
+   *  Deletes all the cells, clears the hashtable,
+   *  resets size, size of cell, number of cells.
+   */
+  virtual void clear();
+
+  /*! Sets the different parameters of the grid
+   *    orig
+   *      The grid origin
+   *    size
+   *      The grid's dimensions
+   *    nb
+   *      The number of cells of the grid
+   */
+  virtual void configure(const Vec3r& orig, const Vec3r& size, unsigned nb);
+  
+  /*! returns the cell whose coordinates
+   *  are pased as argument
+   */
+  virtual Cell* getCell(const Vec3u& p) {
+    Cell* found_cell = NULL;
+    
+    GridHashTable::const_iterator found = _cells.find(p);
+    if (found != _cells.end())
+      found_cell = (*found).second;
+    return found_cell;
+  }
+  
+  /*! Fills the case p with the cell iCell */
+  virtual void fillCell(const Vec3u& p, Cell& cell) {
+    _cells[p] = &cell;
+  }
+
+protected:
+
+  GridHashTable _cells;
+};
+
+#endif // HASHGRID_H
diff --git a/source/blender/freestyle/intern/geometry/Noise.cpp b/source/blender/freestyle/intern/geometry/Noise.cpp
new file mode 100755
index 00000000000..396fc3bbb47
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/Noise.cpp
@@ -0,0 +1,264 @@
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#include "Noise.h"
+# include <stdlib.h>
+# include <stdio.h>
+# include <math.h>
+#include <time.h>
+
+#define MINX		-1000000
+#define MINY		MINX
+#define MINZ		MINX
+#define SCURVE(a) ((a)*(a)*(3.0-2.0*(a)))
+#define REALSCALE ( 2.0 / 65536.0 )
+#define NREALSCALE ( 2.0 / 4096.0 )
+#define HASH3D(a,b,c) hashTable[hashTable[hashTable[(a) & 0xfff] ^ ((b) & 0xfff)] ^ ((c) & 0xfff)]
+#define HASH(a,b,c) (xtab[(xtab[(xtab[(a) & 0xff] ^ (b)) & 0xff] ^ (c)) & 0xff] & 0xff)
+#define INCRSUM(m,s,x,y,z)	((s)*(RTable[m]*0.5		\
+					+ RTable[m+1]*(x)	\
+					+ RTable[m+2]*(y)	\
+					+ RTable[m+3]*(z)))
+#define MAXSIZE		500
+#define nrand()		((float)rand()/(float)RAND_MAX)
+#define seednrand(x)	srand(x*RAND_MAX)
+
+#define BM 0xff
+
+#define N 0x1000
+#define NP 12   /* 2^N */
+#define NM 0xfff
+
+#define s_curve(t) ( t * t * (3. - 2. * t) )
+
+#define lerp(t, a, b) ( a + t * (b - a) )
+
+#define setup(i,b0,b1,r0,r1)\
+	t = i + N;\
+	b0 = ((int)t) & BM;\
+	b1 = (b0+1) & BM;\
+	r0 = t - (int)t;\
+	r1 = r0 - 1.;
+
+void normalize2(float v[2])
+{
+  float s;
+
+  s = sqrt(v[0] * v[0] + v[1] * v[1]);
+  v[0] = v[0] / s;
+  v[1] = v[1] / s;
+}
+
+void normalize3(float v[3])
+{
+  float s;
+
+  s = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
+  v[0] = v[0] / s;
+  v[1] = v[1] / s;
+  v[2] = v[2] / s;
+}
+
+float Noise::turbulence1(float arg, float freq, float amp, unsigned oct)
+{
+  float t;
+  float vec;
+  
+  for (t = 0; oct > 0 && freq > 0; freq *= 2, amp /= 2, --oct) 
+    {
+      vec = freq * arg;
+      t += smoothNoise1(vec) * amp;
+    }
+  return t;
+}
+
+float Noise::turbulence2(Vec2f& v, float freq, float amp, unsigned oct)
+{
+  float t;
+  Vec2f vec;
+  
+  for (t = 0; oct > 0 && freq > 0; freq *= 2, amp /= 2, --oct) 
+    {
+      vec.x() = freq * v.x();
+      vec.y() = freq * v.y();
+      t += smoothNoise2(vec) * amp;
+    }
+  return t;
+}
+
+float Noise::turbulence3(Vec3f& v, float freq, float amp, unsigned oct)
+{
+  float t;
+  Vec3f vec;
+  
+  for (t = 0; oct > 0 && freq > 0; freq *= 2, amp /= 2, --oct) 
+    {
+      vec.x() = freq * v.x();
+      vec.y() = freq * v.y();
+      vec.z() = freq * v.z();
+      t += smoothNoise3(vec) * amp;
+    }
+  return t;
+}
+
+// Noise functions over 1, 2, and 3 dimensions
+
+float Noise::smoothNoise1(float arg)
+{
+  int bx0, bx1;
+  float rx0, rx1, sx, t, u, v, vec;
+
+  vec = arg;
+  setup(vec, bx0,bx1, rx0,rx1);
+
+  sx = s_curve(rx0);
+
+  u = rx0 * g1[ p[ bx0 ] ];
+  v = rx1 * g1[ p[ bx1 ] ];
+
+  return lerp(sx, u, v);
+}
+
+float Noise::smoothNoise2(Vec2f& vec)
+{
+  int bx0, bx1, by0, by1, b00, b10, b01, b11;
+  float rx0, rx1, ry0, ry1, *q, sx, sy, a, b, t, u, v;
+  register int i, j;
+
+  setup(vec.x(), bx0,bx1, rx0,rx1);
+  setup(vec.y(), by0,by1, ry0,ry1);
+
+  i = p[ bx0 ];
+  j = p[ bx1 ];
+
+  b00 = p[ i + by0 ];
+  b10 = p[ j + by0 ];
+  b01 = p[ i + by1 ];
+  b11 = p[ j + by1 ];
+
+  sx = s_curve(rx0);
+  sy = s_curve(ry0);
+
+#define at2(rx,ry) ( rx * q[0] + ry * q[1] )
+
+  q = g2[ b00 ] ; u = at2(rx0,ry0);
+  q = g2[ b10 ] ; v = at2(rx1,ry0);
+  a = lerp(sx, u, v);
+
+  q = g2[ b01 ] ; u = at2(rx0,ry1);
+  q = g2[ b11 ] ; v = at2(rx1,ry1);
+  b = lerp(sx, u, v);
+
+  return lerp(sy, a, b);
+}
+
+float Noise::smoothNoise3(Vec3f& vec)
+{
+  int bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11;
+  float rx0, rx1, ry0, ry1, rz0, rz1, *q, sy, sz, a, b, c, d, t, u, v;
+  register int i, j;
+
+  setup(vec.x(), bx0,bx1, rx0,rx1);
+  setup(vec.y(), by0,by1, ry0,ry1);
+  setup(vec.z(), bz0,bz1, rz0,rz1);
+
+  i = p[ bx0 ];
+  j = p[ bx1 ];
+
+  b00 = p[ i + by0 ];
+  b10 = p[ j + by0 ];
+  b01 = p[ i + by1 ];
+  b11 = p[ j + by1 ];
+
+  t  = s_curve(rx0);
+  sy = s_curve(ry0);
+  sz = s_curve(rz0);
+
+#define at3(rx,ry,rz) ( rx * q[0] + ry * q[1] + rz * q[2] )
+
+  q = g3[ b00 + bz0 ] ;
+  u = at3(rx0,ry0,rz0);
+  q = g3[ b10 + bz0 ] ;
+  v = at3(rx1,ry0,rz0);
+  a = lerp(t, u, v);
+
+  q = g3[ b01 + bz0 ] ;
+  u = at3(rx0,ry1,rz0);
+  q = g3[ b11 + bz0 ] ;
+  v = at3(rx1,ry1,rz0);
+  b = lerp(t, u, v);
+
+  c = lerp(sy, a, b);
+
+  q = g3[ b00 + bz1 ] ;
+  u = at3(rx0,ry0,rz1);
+  q = g3[ b10 + bz1 ] ;
+  v = at3(rx1,ry0,rz1);
+  a = lerp(t, u, v);
+
+  q = g3[ b01 + bz1 ] ;
+  u = at3(rx0,ry1,rz1);
+  q = g3[ b11 + bz1 ] ;
+  v = at3(rx1,ry1,rz1);
+  b = lerp(t, u, v);
+
+  d = lerp(sy, a, b);
+
+  return lerp(sz, c, d);
+}
+
+Noise::Noise()
+{
+  int i, j, k;
+  
+  seednrand(time(NULL));
+  for (i = 0 ; i < _Noise_B_ ; i++) 
+    {
+      p[i] = i;
+
+      g1[i] = (float)((rand() % (_Noise_B_ + _Noise_B_)) - _Noise_B_) / _Noise_B_;
+
+      for (j = 0 ; j < 2 ; j++)
+	g2[i][j] = (float)((rand() % (_Noise_B_ + _Noise_B_)) - _Noise_B_) / _Noise_B_;
+      normalize2(g2[i]);
+
+      for (j = 0 ; j < 3 ; j++)
+	g3[i][j] = (float)((rand() % (_Noise_B_ + _Noise_B_)) - _Noise_B_) / _Noise_B_;
+      normalize3(g3[i]);
+    }
+
+  while (--i) 
+    {
+      k = p[i];
+      p[i] = p[j = rand() % _Noise_B_];
+      p[j] = k;
+    }
+
+  for (i = 0 ; i < _Noise_B_ + 2 ; i++) 
+    {
+      p[_Noise_B_ + i] = p[i];
+      g1[_Noise_B_ + i] = g1[i];
+      for (j = 0 ; j < 2 ; j++)
+	g2[_Noise_B_ + i][j] = g2[i][j];
+      for (j = 0 ; j < 3 ; j++)
+	g3[_Noise_B_ + i][j] = g3[i][j];
+    }
+}
diff --git a/source/blender/freestyle/intern/geometry/Noise.h b/source/blender/freestyle/intern/geometry/Noise.h
new file mode 100755
index 00000000000..00cebbb451e
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/Noise.h
@@ -0,0 +1,77 @@
+//
+//  Filename         : Noise.h
+//  Author(s)        : Emmanuel Turquin
+//  Purpose          : Class to define Perlin noise
+//  Date of creation : 12/01/2004
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef NOISE_H
+# define NOISE_H
+
+
+# include "../system/FreestyleConfig.h"
+# include "Geom.h"
+
+#define _Noise_B_ 0x100
+
+using namespace Geometry;
+using namespace std;
+
+/*! Class to provide Perlin Noise functionalities */
+class LIB_GEOMETRY_EXPORT Noise
+{
+ public:
+
+  /*! Builds a Noise object */
+  Noise();
+  /*! Destructor */
+  ~Noise() {}
+
+  /*! Returns a noise value for a 1D element */
+  float turbulence1(float arg, float freq, float amp, unsigned oct = 4);
+
+  /*! Returns a noise value for a 2D element */
+  float turbulence2(Vec2f& v, float freq, float amp, unsigned oct = 4);
+
+  /*! Returns a noise value for a 3D element */
+  float turbulence3(Vec3f& v, float freq, float amp, unsigned oct = 4);
+
+  /*! Returns a smooth noise value for a 1D element */
+  float smoothNoise1(float arg);
+  /*! Returns a smooth noise value for a 2D element */
+  float smoothNoise2(Vec2f& vec);
+  /*! Returns a smooth noise value for a 3D element */
+  float smoothNoise3(Vec3f& vec);
+
+ private:
+
+  int p[ _Noise_B_ + _Noise_B_ + 2];
+  float g3[ _Noise_B_ + _Noise_B_ + 2][3];
+  float g2[ _Noise_B_ + _Noise_B_ + 2][2];
+  float g1[ _Noise_B_ + _Noise_B_ + 2];
+  int start;
+};
+
+#endif // NOISE_H
diff --git a/source/blender/freestyle/intern/geometry/Polygon.h b/source/blender/freestyle/intern/geometry/Polygon.h
new file mode 100755
index 00000000000..f9c4c78d424
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/Polygon.h
@@ -0,0 +1,215 @@
+//
+//  Filename         : Polygon.h
+//  Author(s)        : Stephane Grabli
+//  Purpose          : Class to define a polygon
+//  Date of creation : 30/07/2002
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef  POLYGON_H
+# define POLYGON_H
+
+# include <vector>
+# include "Geom.h"
+# include "GeomUtils.h"
+
+using namespace std;
+
+namespace Geometry {
+
+template <class Point>
+class Polygon
+{
+ public:
+
+  inline Polygon() {
+    _id = 0;
+    userdata = 0;
+    userdata2 = 0;
+  }
+
+  inline Polygon(const vector<Point>& vertices) {
+    _vertices = vertices;
+    computeBBox();
+    _id = 0;
+    userdata = 0;
+    userdata2 = 0;
+  }
+
+  inline Polygon(const Polygon<Point>& poly) {
+    Point p;
+    for(typename vector<Point>::const_iterator it = poly.getVertices().begin();
+	it != poly.getVertices().end();
+	it++) {
+      p = *it;
+      _vertices.push_back(p);
+    }
+
+    _id = poly.getId();
+    poly.getBBox(_min, _max);
+    userdata = 0;
+    userdata2 = 0;
+  }
+  
+  virtual ~Polygon() {}
+
+  //
+  // Accessors
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  inline const vector<Point>& getVertices() const {
+    return _vertices;
+  }
+
+  inline void getBBox(Point& min, Point& max) const {
+    min = _min;
+    max = _max;
+  }
+
+  inline Point& getBBoxCenter()
+  {
+    Point result;
+    result = (_min + _max) / 2;
+    return result;
+  }
+
+  inline Point& getCenter() {
+    Point result;
+    for (typename vector<Point>::iterator it = _vertices.begin();
+	 it != _vertices.end();
+	 it++)
+      result += *it;
+    result /= _vertices.size();
+    return result;
+  }
+
+  inline unsigned getId() const {
+    return _id;
+  }
+
+  //
+  // Modifiers
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  inline void setVertices(const vector<Point>& vertices) {
+    _vertices.clear();
+    Point p;
+    for (typename vector<Point>::const_iterator it = vertices.begin();
+	 it != vertices.end();
+	 it++) {
+      p = *it;
+      _vertices.push_back(p);
+    }
+    computeBBox();
+  }
+
+  inline void setId(unsigned id) {
+    _id = id;
+  }
+
+  //
+  // Other methods
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  inline void computeBBox() {
+    if(_vertices.empty())
+      return;
+    
+    _max = _vertices[0];
+    _min = _vertices[0];
+
+    for(typename vector<Point>::iterator it = _vertices.begin();
+	it != _vertices.end();
+        it++) {
+      for(unsigned i = 0; i < Point::dim(); i++) {
+	if((*it)[i] > _max[i])
+	  _max[i] = (*it)[i];
+	if((*it)[i] < _min[i])
+	  _min[i] = (*it)[i];
+      }
+    }
+  }
+  
+  // FIXME Is it possible to get rid of userdatas ?
+  void* userdata;
+  void* userdata2; // Used during ray casting
+  
+ protected:
+  
+  vector<Point>	_vertices;
+  Point		_min;
+  Point		_max;
+  unsigned	_id;
+};
+
+
+//
+// Polygon3r class
+//
+///////////////////////////////////////////////////////////////////////////////
+
+class Polygon3r : public Polygon<Vec3r>
+{
+ public:
+
+  inline Polygon3r() : Polygon<Vec3r>() {}
+
+  inline Polygon3r(const vector<Vec3r>& vertices,
+		   const Vec3r& normal) : Polygon<Vec3r>(vertices) {
+    setNormal(normal);
+  }
+
+  inline Polygon3r(const Polygon3r& poly) : Polygon<Vec3r>(poly) {}
+  
+  virtual ~Polygon3r() {}
+
+  void setNormal(const Vec3r& normal) {
+    _normal = normal;
+  }
+
+  Vec3r getNormal() const {
+    return _normal;
+  }
+
+  /*! Check whether the Polygon intersects with the ray or not */
+  inline bool rayIntersect(Vec3r& orig, Vec3r& dir,
+			   real& t, real& u, real& v, real epsilon = M_EPSILON) {
+    //    if (_vertices.size() < 3)
+    //      return false;
+    return GeomUtils::intersectRayTriangle(orig, dir,
+					   _vertices[0], _vertices[1], _vertices[2],
+					   t, u, v, epsilon);
+  }
+
+ private:
+
+  Vec3r	_normal;
+};
+
+} // end of namespace Geometry
+
+#endif // POLYGON_H
diff --git a/source/blender/freestyle/intern/geometry/SweepLine.h b/source/blender/freestyle/intern/geometry/SweepLine.h
new file mode 100755
index 00000000000..e3fb4ad8c0c
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/SweepLine.h
@@ -0,0 +1,334 @@
+//
+//  Filename         : SweepLine.h
+//  Author(s)        : Stephane Grabli
+//  Purpose          : Class to define a Sweep Line
+//  Date of creation : 29/08/2002
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef  SWEEPLINE_H
+# define SWEEPLINE_H
+
+# include <list>
+# include <vector>
+
+/*! Class to define the intersection berween two segments*/
+template<class Edge>
+class Intersection
+{
+public:
+
+  template<class EdgeClass>
+  Intersection(EdgeClass* eA, real ta, EdgeClass* eB, real tb)
+  {
+    EdgeA = eA;
+    EdgeB = eB;
+    tA = ta;
+    tB = tb;
+    userdata = 0;
+  }
+
+  Intersection(const Intersection& iBrother)
+  {
+    EdgeA = iBrother.EdgeA;
+    EdgeB = iBrother.EdgeB;
+    tA = iBrother.tA;
+    tB = iBrother.tB;
+    userdata = 0;
+  }
+
+  /*! returns the parameter giving the 
+   *  intersection, for the edge iEdge 
+   */
+  real getParameter(Edge *iEdge) 
+  {
+    if(iEdge == EdgeA)
+      return tA;
+    if(iEdge == EdgeB)
+      return tB;
+    return 0;
+  }
+
+public:
+  void * userdata; // FIXME
+
+  Edge *EdgeA; // first segment
+  Edge *EdgeB; // second segment
+  real tA; // parameter defining the intersection point with respect to the segment EdgeA.
+  real tB; // parameter defining the intersection point with respect to the segment EdgeB.
+};
+
+
+
+
+
+
+
+
+template<class T, class Point>
+class Segment
+{
+  public:
+    Segment()
+    {
+    }
+    Segment(T& s, const Point& iA, const Point& iB) 
+    {
+      _edge = s;
+      if(iA < iB)
+      {
+        A = iA;
+        B = iB;
+        _order = true;
+      }
+      else
+      {
+        A = iB;
+        B = iA;
+        _order = false;
+      }
+    }
+
+    Segment(Segment<T,Point>& iBrother)
+    {
+      _edge = iBrother.edge();
+      A = iBrother.A;
+      B = iBrother.B;
+      _Intersections = iBrother._Intersections;
+      _order = iBrother._order;
+    }
+
+    Segment(const Segment<T,Point>& iBrother)
+    {
+      _edge = iBrother._edge;
+      A = iBrother.A;
+      B = iBrother.B;
+      _Intersections = iBrother._Intersections;
+      _order = iBrother._order;
+    }
+
+    ~Segment() {
+      _Intersections.clear();
+    }
+
+    inline Point operator[](const unsigned short int& i) const
+    {
+      return i%2==0 ? A : B;
+    }
+
+    inline bool operator==(const Segment<T,Point>& iBrother)
+    {
+      if(_edge == iBrother._edge)
+        return true;
+
+      return false;
+    }
+  
+    /* Adds an intersection for this segment */
+    inline void AddIntersection(Intersection<Segment<T,Point> > *i) {_Intersections.push_back(i);}
+
+    /*! Checks for a common vertex with another edge */
+    inline bool CommonVertex(const Segment<T,Point>& S, Point& CP)
+    {
+      if((A == S[0]) || (A == S[1]))
+      {
+        CP = A;
+        return true;
+      }
+      if((B == S[0]) || (B == S[1]))
+      {
+        CP = B;
+        return true;
+      }
+
+      return false;
+    }
+
+    inline vector<Intersection<Segment<T,Point> >*>& intersections() {return _Intersections;}
+    inline bool order() {return _order;}
+    inline T& edge() {return _edge;}
+
+  private:
+    T _edge;
+    Point A;
+    Point B;
+    std::vector<Intersection<Segment<T,Point> >*> _Intersections; // list of intersections parameters
+    bool _order; // true if A and B are in the same order than _edge.A and _edge.B. false otherwise.
+};
+
+/*! defines a binary function that can be overload 
+ *  by the user to specify at each condition 
+ *  the intersection between 2 edges must be computed
+ */
+template<class T1, class T2>
+struct binary_rule 
+{
+  binary_rule() {}
+  template<class T3,class T4>
+  binary_rule(const binary_rule<T3,T4>& brother) {}
+  virtual ~binary_rule() {}
+
+  virtual bool operator()(T1&, T2&)
+  {
+    return true;
+  }
+};
+
+
+template<class T,class Point>
+class SweepLine 
+{
+public:
+
+  SweepLine() {}
+  ~SweepLine() 
+  {
+    for(typename vector<Intersection<Segment<T,Point> >*>::iterator i=_Intersections.begin(),iend=_Intersections.end();
+    i!=iend;
+    i++)
+    {
+      delete (*i);
+    }
+    _Intersections.clear();
+
+    for(typename vector<Segment<T,Point>* >::iterator ie=_IntersectedEdges.begin(),ieend=_IntersectedEdges.end();
+    ie!=ieend;
+    ie++)
+    {
+      delete (*ie);
+    }
+    _IntersectedEdges.clear();
+
+    _set.clear();
+  }
+
+  
+  inline void process(Point& p, 
+                      vector<Segment<T,Point>*>& segments, 
+                      binary_rule<Segment<T,Point>,Segment<T,Point> >& binrule
+                      //binary_rule<Segment<T,Point>,Segment<T,Point> >& binrule = binary_rule<Segment<T,Point>,Segment<T,Point> >()
+                      ) 
+  {
+    // first we remove the segments that need to be removed and then 
+    // we add the segments to add
+    vector<Segment<T,Point>*> toadd;
+    typename vector<Segment<T,Point>*>::iterator s, send;
+    for(s=segments.begin(), send=segments.end();
+    s!=send;
+    s++)
+    {
+       if(p == (*(*s))[0])
+        toadd.push_back((*s));
+      else
+        remove((*s));
+    }
+    for(s=toadd.begin(), send=toadd.end();
+    s!=send;
+    s++)
+    {
+      add((*s), binrule);
+    }
+  }
+  
+  inline void add(Segment<T,Point>* S,
+                  binary_rule<Segment<T,Point>,Segment<T,Point> >& binrule
+                  //binary_rule<Segment<T,Point>,Segment<T,Point> >& binrule = binary_rule<Segment<T,Point>, Segment<T,Point> >()
+                  )
+  {
+    real t,u;
+    Point CP;
+    Vec2r v0, v1, v2, v3;
+    if(true == S->order())
+    {
+      v0[0] = ((*S)[0])[0];
+      v0[1] = ((*S)[0])[1];
+      v1[0] = ((*S)[1])[0];
+      v1[1] = ((*S)[1])[1];
+    }
+    else
+    {
+      v1[0] = ((*S)[0])[0];
+      v1[1] = ((*S)[0])[1];
+      v0[0] = ((*S)[1])[0];
+      v0[1] = ((*S)[1])[1];
+    }
+    for(typename std::list<Segment<T,Point>* >::iterator s=_set.begin(), send=_set.end();
+    s!=send;
+    s++)
+    {
+      Segment<T,Point>* currentS = (*s);
+      if(true != binrule(*S, *currentS))
+        continue;
+
+      if(true == currentS->order())
+      {
+        v2[0] = ((*currentS)[0])[0];
+	v2[1] = ((*currentS)[0])[1];
+        v3[0] = ((*currentS)[1])[0];
+	v3[1] = ((*currentS)[1])[1];
+      }
+      else
+      {
+        v3[0] = ((*currentS)[0])[0];
+	v3[1] = ((*currentS)[0])[1];
+        v2[0] = ((*currentS)[1])[0];
+	v2[1] = ((*currentS)[1])[1];
+      }
+      if(S->CommonVertex(*currentS, CP))
+        continue; // the two edges have a common vertex->no need to check
+      
+      if(GeomUtils::intersect2dSeg2dSegParametric(v0, v1, v2, v3, t, u))
+      {
+        // create the intersection
+        Intersection<Segment<T,Point> > * inter = new Intersection<Segment<T,Point> >(S,t,currentS,u);
+        // add it to the intersections list
+        _Intersections.push_back(inter);
+        // add this intersection to the first edge intersections list
+        S->AddIntersection(inter);
+        // add this intersection to the second edge intersections list
+        currentS->AddIntersection(inter);        
+      }
+    }
+    // add the added segment to the list of active segments
+    _set.push_back(S);
+  }
+  
+  inline void remove(Segment<T,Point>* s) 
+  {
+    if(s->intersections().size() > 0)
+      _IntersectedEdges.push_back(s);
+    _set.remove(s);
+  }
+
+  vector<Segment<T,Point>* >& intersectedEdges() {return _IntersectedEdges;}
+  vector<Intersection<Segment<T,Point> >*>& intersections() {return _Intersections;}
+  
+
+private:
+  std::list<Segment<T,Point>* > _set; // set of active edges for a given position of the sweep line
+  std::vector<Segment<T,Point>* > _IntersectedEdges; // the list of intersected edges
+  std::vector<Intersection<Segment<T,Point> >*> _Intersections; // the list of all intersections.
+};
+
+#endif // SWEEPLINE_H
diff --git a/source/blender/freestyle/intern/geometry/VecMat.h b/source/blender/freestyle/intern/geometry/VecMat.h
new file mode 100755
index 00000000000..9bbec3b1349
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/VecMat.h
@@ -0,0 +1,899 @@
+//
+//  Filename         : VecMat.h
+//  Author(s)        : Sylvain Paris
+//                     Emmanuel Turquin
+//                     Stephane Grabli 
+//  Purpose          : Vectors and Matrices definition and manipulation
+//  Date of creation : 12/06/2003
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef  VECMAT_H
+# define VECMAT_H
+
+# include <math.h>
+# include <vector>
+# include <iostream>
+
+namespace VecMat {
+
+  namespace Internal {
+
+    template <bool B>
+    struct is_false {};
+
+    template <>
+    struct is_false<false> {
+      static inline void ensure() {}
+    };
+
+  } // end of namespace Internal
+
+  //
+  //  Vector class
+  //    - T: value type
+  //    - N: dimension
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  template <class T, unsigned N>
+  class Vec
+  {
+  public:
+
+    typedef T value_type;
+    
+    // constructors
+
+    inline Vec() {
+      for (unsigned i = 0; i < N; i++)
+	this->_coord[i] = 0;
+    }
+
+    ~Vec() {
+      Internal::is_false<(N == 0)>::ensure();
+    }
+
+    template <class U>
+    explicit inline Vec(const U tab[N]) {
+      for (unsigned i = 0; i < N; i++)
+	this->_coord[i] = (T)tab[i];
+    }
+
+    template <class U>
+    explicit inline Vec(const std::vector<U>& tab) {
+      for (unsigned i = 0; i < N; i++)
+	this->_coord[i] = (T)tab[i];
+    }
+
+    template <class U>
+    explicit inline Vec(const Vec<U, N>& v) {
+      for (unsigned i = 0; i < N; i++)
+        this->_coord[i] = (T)v[i];
+    }
+
+    // accessors
+
+    inline value_type  operator[](const unsigned i) const {
+      return this->_coord[i];
+    }
+
+    inline value_type& operator[](const unsigned i) {
+      return this->_coord[i];
+    }
+
+    static inline unsigned dim() {
+      return N;
+    }
+    
+    // various useful methods
+
+    inline value_type norm() const {
+      return (T)sqrt((float)squareNorm());
+    }
+
+    inline value_type squareNorm() const {
+      return (*this) * (*this);
+    }
+
+    inline Vec<T, N>& normalize() {
+      value_type n = norm();
+      for (unsigned i = 0; i < N; i++)
+	this->_coord[i] /= n;
+      return *this;
+    }
+
+    inline Vec<T, N>& normalizeSafe() {
+      value_type n = norm();
+      if (n)
+	for (unsigned i=0; i < N; i++)
+	  this->_coord[i] /= n;
+      return *this;
+    }
+
+    // classical operators
+    inline Vec<T, N> operator+(const Vec<T, N>& v) const{
+    Vec<T, N> res(v);
+    res += *this;
+    return res;
+    } 
+
+    inline Vec<T, N> operator-(const Vec<T,N>& v) const{
+    Vec<T, N> res(*this);
+    res -= v;
+    return res;
+    } 
+
+    inline Vec<T, N> operator*(const typename Vec<T,N>::value_type r) const{
+    Vec<T, N> res(*this);
+    res *= r;
+    return res;
+    } 
+
+    inline Vec<T, N> operator/(const typename Vec<T,N>::value_type r) const{
+    Vec<T, N> res(*this);
+    if (r)
+      res /= r;
+    return res;
+    } 
+
+    // dot product
+    inline value_type operator*(const Vec<T, N>& v) const{
+    value_type sum = 0;
+    for (unsigned i = 0; i < N; i++)
+      sum += (*this)[i] * v[i];
+    return sum;
+    } 
+    
+    template <class U>
+    inline Vec<T, N>& operator=(const Vec<U, N>& v) {
+      if (this != &v)
+	for (unsigned i = 0; i < N; i++)
+	  this->_coord[i] = (T)v[i];
+      return *this;
+    }
+
+    template <class U>
+    inline Vec<T, N>& operator+=(const Vec<U, N>& v) {
+      for (unsigned i = 0 ; i < N; i++)
+	this->_coord[i] += (T)v[i];
+      return *this;
+    }
+
+    template <class U>
+    inline Vec<T, N>& operator-=(const Vec<U, N>& v) {
+      for (unsigned i = 0 ; i < N; i++)
+	this->_coord[i] -= (T)v[i];
+      return *this;
+    }
+
+    template <class U>
+    inline Vec<T, N>& operator*=(const U r) {
+      for (unsigned i = 0 ; i < N; i++)
+	this->_coord[i] *= r;
+      return *this;
+    }
+
+    template <class U>
+    inline Vec<T, N>& operator/=(const U r) {
+      if (r)
+	for (unsigned i = 0 ; i < N; i++)
+	  this->_coord[i] /= r;
+      return *this;
+    }
+
+
+    inline bool operator==(const Vec<T, N>& v) const {
+      for(unsigned i = 0; i < N; i++)
+	if (this->_coord[i] != v[i])
+	  return false;
+      return true;
+    }
+
+    inline bool operator!=(const Vec<T, N>& v) const {
+      for(unsigned i = 0; i < N; i++)
+	if (this->_coord[i] != v[i])
+	  return true;
+      return false;
+    }
+
+    inline bool operator<(const Vec<T, N>& v) const {
+      for (unsigned i = 0; i<N; i++) {
+	if (this->_coord[i] < v[i])
+	  return true;
+	if (this->_coord[i] > v[i])
+	  return false;
+	if (this->_coord[i] == v[i])
+	  continue;
+      }
+      return false;  
+    }
+
+    inline bool operator>(const Vec<T, N>& v) const {
+      for (unsigned i=0; i<N; i++) {
+	if(this->_coord[i] > v[i])
+	  return true;
+	if(this->_coord[i] < v[i])
+	  return false;
+	if(this->_coord[i] == v[i])
+	  continue;
+      }
+      return false;  
+    }
+
+  protected:
+
+    value_type _coord[N];
+    enum {
+      _dim = N,
+    };
+  };
+
+
+  //
+  //  Vec2 class (2D Vector)
+  //    - T: value type
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  template <class T>
+  class Vec2 : public Vec<T, 2>
+  {
+  public:
+
+    typedef typename Vec<T, 2>::value_type value_type;
+
+    inline Vec2() : Vec<T, 2>() {}
+
+    template <class U>
+    explicit inline Vec2(const U tab[2]) : Vec<T, 2>(tab) {}
+
+    template <class U>
+    explicit inline Vec2(const std::vector<U>& tab) : Vec<T, 2>(tab) {}
+
+    template <class U>
+    inline Vec2(const Vec<U, 2>& v) : Vec<T, 2>(v) {}
+
+    inline Vec2(const value_type x,
+		const value_type y = 0) : Vec<T, 2>() {
+      this->_coord[0] = (T)x;
+      this->_coord[1] = (T)y;
+    }
+    
+    inline value_type x() const {
+      return this->_coord[0];
+    }
+
+    inline value_type& x() {
+      return this->_coord[0];
+    }
+
+    inline value_type y() const {
+      return this->_coord[1];
+    }
+
+    inline value_type& y() {
+      return this->_coord[1];
+    }
+
+    inline void setX(const value_type v) {
+      this->_coord[0] = v;
+    }
+
+    inline void setY(const value_type v) {
+      this->_coord[1] = v;
+    }
+
+    // FIXME: hack swig -- no choice
+    inline Vec2<T> operator+(const Vec2<T>& v) const{
+    Vec2<T> res(v);
+    res += *this;
+    return res;
+    } 
+
+    inline Vec2<T> operator-(const Vec2<T>& v) const{
+    Vec2<T> res(*this);
+    res -= v;
+    return res;
+    } 
+
+    inline Vec2<T> operator*(const value_type r) const{
+    Vec2<T> res(*this);
+    res *= r;
+    return res;
+    } 
+
+    inline Vec2<T> operator/(const value_type r) const{
+    Vec2<T> res(*this);
+    if (r)
+      res /= r;
+    return res;
+    } 
+
+    // dot product
+    inline value_type operator*(const Vec2<T>& v) const{
+    value_type sum = 0;
+    for (unsigned i = 0; i < 2; i++)
+      sum += (*this)[i] * v[i];
+    return sum;
+    } 
+  };
+
+
+  //
+  //  HVec3 class (3D Vector in homogeneous coordinates)
+  //    - T: value type
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  template <class T>
+  class HVec3 : public Vec<T, 4>
+  {
+  public:
+
+    typedef typename Vec<T, 4>::value_type value_type;
+
+    inline HVec3() : Vec<T, 4>() {}
+
+    template <class U>
+    explicit inline HVec3(const U tab[4]) : Vec<T, 4>(tab) {}
+
+    template <class U>
+    explicit inline HVec3(const std::vector<U>& tab) : Vec<T, 4>(tab) {}
+
+    template<class U>
+    inline HVec3(const Vec<U, 4>& v) : Vec<T, 4>(v) {}
+  
+    inline HVec3(const value_type sx,
+		 const value_type sy = 0,
+		 const value_type sz = 0,
+		 const value_type s = 1) {
+      this->_coord[0] = sx;
+      this->_coord[1] = sy;
+      this->_coord[2] = sz;
+      this->_coord[3] = s;
+    }
+    
+    template <class U>
+    inline HVec3(const Vec<U, 3>& sv,
+		 const U s = 1) {
+      this->_coord[0] = (T)sv[0];
+      this->_coord[1] = (T)sv[1];
+      this->_coord[2] = (T)sv[2];
+      this->_coord[3] = (T)s;
+    }
+    
+    inline value_type sx() const {
+      return this->_coord[0];
+    }
+
+    inline value_type& sx() {
+      return this->_coord[0];
+    }
+
+    inline value_type sy() const {
+      return this->_coord[1];
+    }
+
+    inline value_type& sy() {
+      return this->_coord[1];
+    }
+
+    inline value_type sz() const {
+      return this->_coord[2];
+    }
+
+    inline value_type& sz() {
+      return this->_coord[2];
+    }
+
+    inline value_type s() const {
+      return this->_coord[3];
+    }
+
+    inline value_type& s() {
+      return this->_coord[3];
+    }
+
+    // Acces to non-homogeneous coordinates in 3D
+
+    inline value_type x() const {
+      return this->_coord[0] / this->_coord[3];
+    }
+
+    inline value_type y() const {
+      return this->_coord[1] / this->_coord[3];
+    }
+
+    inline value_type z() const {
+      return this->_coord[2] / this->_coord[3];
+    }
+  };
+
+
+  //
+  //  Vec3 class (3D Vec)
+  //    - T: value type
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  template <class T>
+  class Vec3 : public Vec<T, 3>
+  {
+  public:
+    
+    typedef typename Vec<T, 3>::value_type value_type;
+
+    inline Vec3() : Vec<T, 3>() {}
+
+    template <class U>
+    explicit inline Vec3(const U tab[3]) : Vec<T, 3>(tab) {}
+
+    template <class U>
+    explicit inline Vec3(const std::vector<U>& tab) : Vec<T, 3>(tab) {}
+
+    template<class U>
+    inline Vec3(const Vec<U, 3>& v) : Vec<T, 3>(v) {}
+  
+    template<class U>
+      inline Vec3(const HVec3<U>& v) {
+      this->_coord[0] = (T)v.x();
+      this->_coord[1] = (T)v.y();
+      this->_coord[2] = (T)v.z();
+    }
+
+    inline Vec3(const value_type x,
+		const value_type y = 0,
+		const value_type z = 0) : Vec<T, 3>() {
+      this->_coord[0] = x;
+      this->_coord[1] = y;
+      this->_coord[2] = z;
+    }
+    
+    inline value_type x() const {
+      return this->_coord[0];
+    }
+
+    inline value_type& x() {
+      return this->_coord[0];
+    }
+
+    inline value_type y() const {
+      return this->_coord[1];
+    }
+
+    inline value_type& y() {
+      return this->_coord[1];
+    }
+
+    inline value_type z() const {
+      return this->_coord[2];
+    }
+
+    inline value_type& z() {
+      return this->_coord[2];
+    }
+
+    inline void setX(const value_type v) {
+      this->_coord[0] = v;
+    }
+
+    inline void setY(const value_type v) {
+      this->_coord[1] = v;
+    }
+
+    inline void setZ(const value_type v) {
+      this->_coord[2] = v;
+    }
+
+    // classical operators
+    // FIXME: hack swig -- no choice
+    inline Vec3<T> operator+(const Vec3<T>& v) const{
+    Vec3<T> res(v);
+    res += *this;
+    return res;
+    } 
+
+    inline Vec3<T> operator-(const Vec3<T>& v) const{
+    Vec3<T> res(*this);
+    res -= v;
+    return res;
+    } 
+
+    inline Vec3<T> operator*(const value_type r) const{
+    Vec3<T> res(*this);
+    res *= r;
+    return res;
+    } 
+
+    inline Vec3<T> operator/(const value_type r) const{
+    Vec3<T> res(*this);
+    if (r)
+      res /= r;
+    return res;
+    } 
+
+    // dot product
+    inline value_type operator*(const Vec3<T>& v) const{
+    value_type sum = 0;
+    for (unsigned i = 0; i < 3; i++)
+      sum += (*this)[i] * v[i];
+    return sum;
+    } 
+
+    // cross product for 3D Vectors
+    // FIXME: hack swig -- no choice
+    inline Vec3<T> operator^(const Vec3<T>& v) const{
+    Vec3<T> res((*this)[1] * v[2] - (*this)[2] * v[1],
+		(*this)[2] * v[0] - (*this)[0] * v[2],
+		(*this)[0] * v[1] - (*this)[1] * v[0]);
+    return res;
+    } 
+
+    // cross product for 3D Vectors
+    template <typename U>
+    inline Vec3<T> operator^(const Vec<U, 3>& v) const{
+    Vec3<T> res((*this)[1] * v[2] - (*this)[2] * v[1],
+		(*this)[2] * v[0] - (*this)[0] * v[2],
+		(*this)[0] * v[1] - (*this)[1] * v[0]);
+    return res;
+    } 
+  };
+
+
+  //
+  //  Matrix class
+  //    - T: value type
+  //    - M: rows
+  //    - N: cols
+  //
+  /////////////////////////////////////////////////////////////////////////////
+  
+  // Dirty, but icc under Windows needs this
+# define _SIZE (M * N)
+
+  template <class T, unsigned M, unsigned N>
+  class Matrix
+  {
+  public:
+
+    typedef T value_type;
+    
+    inline Matrix() {
+      for (unsigned i = 0; i < _SIZE; i++)
+	this->_coord[i] = 0;
+    }
+
+    ~Matrix() {
+      Internal::is_false<(M == 0)>::ensure();
+      Internal::is_false<(N == 0)>::ensure();
+    }
+
+    template <class U>
+    explicit inline Matrix(const U tab[_SIZE]) {
+      for (unsigned i = 0; i < _SIZE; i++)
+	this->_coord[i] = tab[i];
+    }
+
+    template <class U>
+    explicit inline Matrix(const std::vector<U>& tab) {
+      for (unsigned i = 0; i < _SIZE; i++)
+	this->_coord[i] = tab[i];
+    }
+
+    template <class U>
+    inline Matrix(const Matrix<U, M, N>& m) {
+      for (unsigned i = 0; i < M; i++)
+	for (unsigned j = 0; j < N; j++)
+	  this->_coord[i * N + j] = (T)m(i, j);
+    }
+
+    inline value_type operator()(const unsigned i, const unsigned j) const {
+      return this->_coord[i * N + j];
+    }
+
+    inline value_type& operator()(const unsigned i, const unsigned j) {
+      return this->_coord[i * N + j];
+    }
+
+    static inline unsigned rows() {
+      return M;
+    }
+
+    static inline unsigned cols() {
+      return N;
+    }
+
+    inline Matrix<T, M, N>& transpose() const {
+      Matrix<T, N, M> res;
+      for (unsigned i = 0; i < M; i++)
+	for (unsigned j = 0; j < N; j++)
+	  res(j,i) = this->_coord[i * N + j];
+      return res;
+    }
+
+    template <class U>
+    inline Matrix<T, M, N>& operator=(const Matrix<U, M, N>& m) {
+      if (this != &m)
+	for (unsigned i = 0; i < M; i++)
+	  for (unsigned j = 0; j < N; j++)
+	    this->_coord[i * N + j] = (T)m(i, j);
+      return *this;
+    }
+
+    template <class U>
+    inline Matrix<T, M, N>& operator+=(const Matrix<U, M, N>& m) {
+      for (unsigned i = 0; i < M; i++)
+	for (unsigned j = 0; j < N; j++)
+	  this->_coord[i * N + j] += (T)m(i, j);
+      return *this;
+    }
+
+    template <class U>
+    inline Matrix<T, M, N>& operator-=(const Matrix<U, M, N>& m) {
+      for (unsigned i = 0; i < M; i++)
+	for (unsigned j = 0; j < N; j++)
+	  this->_coord[i * N + j] -= (T)m(i, j);
+      return *this;
+    }
+
+    template <class U>
+    inline Matrix<T, M, N>& operator*=(const U lambda) {
+      for (unsigned i = 0; i < M; i++)
+	for (unsigned j = 0; j < N; j++)
+	  this->_coord[i * N + j] *= lambda;
+      return *this;
+    }
+
+    template <class U>
+    inline Matrix<T, M, N>& operator/=(const U lambda) {
+      if (lambda)
+	for (unsigned i = 0; i < M; i++)
+	  for (unsigned j = 0; j < N; j++)
+	    this->_coord[i * N + j] /= lambda;
+      return *this;
+    }
+
+  protected:
+
+    value_type _coord[_SIZE];
+  };
+
+
+  //
+  //  SquareMatrix class
+  //    - T: value type
+  //    - N: rows & cols
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  // Dirty, but icc under Windows needs this
+# define __SIZE (N * N)
+
+  template <class T, unsigned N>
+  class SquareMatrix : public Matrix<T, N, N>
+  {
+  public:
+
+    typedef T value_type;
+  
+    inline SquareMatrix() : Matrix<T, N, N>() {}
+
+    template <class U>
+    explicit inline SquareMatrix(const U tab[__SIZE]) : Matrix<T, N, N>(tab) {}
+
+    template <class U>
+    explicit inline SquareMatrix(const std::vector<U>& tab) : Matrix<T, N, N>(tab) {}
+
+    template <class U>
+    inline SquareMatrix(const Matrix<U, N, N>& m) : Matrix<T, N, N>(m) {}
+
+    static inline SquareMatrix<T, N> identity() {
+      SquareMatrix<T, N> res;
+      for (unsigned i = 0; i < N; i++)
+	res(i, i) = 1;
+      return res;
+    }
+  };
+
+
+  //
+  // Vector external functions
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  //  template <class T, unsigned N>
+  //  inline Vec<T, N> operator+(const Vec<T, N>& v1,
+  //				const Vec<T, N>& v2) {
+  //    Vec<T, N> res(v1);
+  //    res += v2;
+  //    return res;
+  //  }
+  //
+  //  template <class T, unsigned N>
+  //  inline Vec<T, N> operator-(const Vec<T, N>& v1,
+  //				const Vec<T, N>& v2) {
+  //    Vec<T, N> res(v1);
+  //    res -= v2;
+  //    return res;
+  //  }
+  //  template <class T, unsigned N>
+  //    inline Vec<T, N> operator*(const Vec<T, N>& v,
+  //				const typename Vec<T, N>::value_type r) {
+  //    Vec<T, N> res(v);
+  //    res *= r;
+  //    return res;
+  //  }
+  
+  template <class T, unsigned N>
+    inline Vec<T, N> operator*(const typename Vec<T, N>::value_type r,
+				const Vec<T, N>& v) {
+    Vec<T, N> res(v);
+    res *= r;
+    return res;
+  }
+  //
+  //  template <class T, unsigned N>
+  //  inline Vec<T, N> operator/(const Vec<T, N>& v,
+  //				const typename Vec<T, N>::value_type r) {
+  //    Vec<T, N> res(v);
+  //    if (r)
+  //      res /= r;
+  //    return res;
+  //  }
+  //
+  // dot product
+  //  template <class T, unsigned N>
+  //    inline typename Vec<T, N>::value_type operator*(const Vec<T, N>& v1,
+  //    const Vec<T, N>& v2) {
+  //    typename Vec<T, N>::value_type sum = 0;
+  //    for (unsigned i = 0; i < N; i++)
+  //      sum += v1[i] * v2[i];
+  //    return sum;
+  //  }
+  //
+  //  // cross product for 3D Vectors
+  //  template <typename T>
+  //  inline Vec3<T> operator^(const Vec<T, 3>& v1,
+  //			   const Vec<T, 3>& v2) {
+  //    Vec3<T> res(v1[1] * v2[2] - v1[2] * v2[1],
+  //		v1[2] * v2[0] - v1[0] * v2[2],
+  //		v1[0] * v2[1] - v1[1] * v2[0]);
+  //    return res;
+  //  }
+
+  // stream operator
+  template <class T, unsigned N>
+  inline std::ostream& operator<<(std::ostream& s,
+				  const Vec<T, N>& v) {
+    unsigned i;
+    s << "[";
+    for (i = 0; i < N - 1; i++)
+      s << v[i] << ", ";
+    s << v[i] << "]";
+    return s;
+  }
+
+
+  //
+  // Matrix external functions
+  //
+  /////////////////////////////////////////////////////////////////////////////
+
+  template <class T, unsigned M, unsigned N>
+  inline Matrix<T, M, N>
+  operator+(const Matrix<T, M, N>& m1,
+	    const Matrix<T, M, N>& m2) {
+    Matrix<T, M, N> res(m1);
+    res += m2;
+    return res;
+  }
+
+  template <class T, unsigned M, unsigned N>
+  inline Matrix<T, M, N>
+  operator-(const Matrix<T, M, N>& m1,
+	    const Matrix<T, M, N>& m2) {
+    Matrix<T, M, N> res(m1);
+    res -= m2;
+    return res;
+  }
+
+  template <class T, unsigned M, unsigned N>
+  inline Matrix<T, M, N>
+  operator*(const Matrix<T, M, N>& m1,
+	    const typename Matrix<T, M, N>::value_type lambda) {
+    Matrix<T, M, N> res(m1);
+    res *= lambda;
+    return res;
+  }
+
+  template <class T, unsigned M, unsigned N>
+  inline Matrix<T, M, N>
+  operator*(const typename Matrix<T, M, N>::value_type lambda,
+	    const Matrix<T, M, N>& m1) {
+    Matrix<T, M, N> res(m1);
+    res *= lambda;
+    return res;
+  }
+
+  template <class T, unsigned M, unsigned N>
+  inline Matrix<T, M, N>
+  operator/(const Matrix<T, M, N>& m1,
+	    const typename Matrix<T, M, N>::value_type lambda) {
+    Matrix<T, M, N> res(m1);
+    res /= lambda;
+    return res;
+  }
+
+  template <class T, unsigned M, unsigned N, unsigned P>
+  inline Matrix<T, M, P>
+  operator*(const Matrix<T, M, N>& m1,
+	    const Matrix<T, N, P>& m2) {
+    unsigned i, j, k;
+    Matrix<T, M, P> res;
+    typename  Matrix<T, N, P>::value_type scale;
+  
+    for (j = 0; j < P; j++) {
+      for (k = 0; k < N; k++) {
+	scale = m2(k, j);
+	for (i = 0; i < N; i++)
+	  res(i, j) += m1(i, k) * scale;
+      }
+    }
+    return res;
+  }
+
+  template <class T, unsigned M, unsigned N>
+  inline Vec<T, M>
+  operator*(const Matrix<T, M, N>& m,
+	    const Vec<T, N>& v) {
+  
+    Vec<T, M> res;
+    typename Matrix<T, M, N>::value_type scale;
+  
+    for (unsigned j = 0; j < M; j++) {
+      scale = v[j];
+      for (unsigned i = 0; i < N; i++)
+	res[i] += m(i, j) * scale;
+    }
+    return res;
+  }
+
+  // stream operator
+  template <class T, unsigned M, unsigned N>
+  inline std::ostream& operator<<(std::ostream& s,
+				  const Matrix<T, M, N>& m) {
+    unsigned i, j;
+    for (i = 0; i < M; i++) {
+      s << "[";
+      for (j = 0; j < N - 1; j++)
+	s << m(i, j) << ", ";
+      s << m(i, j) << "]" << std::endl;
+    }
+    return s;
+  }
+
+} // end of namespace VecMat
+
+#endif // VECMAT_H
diff --git a/source/blender/freestyle/intern/geometry/geometry.pro b/source/blender/freestyle/intern/geometry/geometry.pro
new file mode 100755
index 00000000000..a63aa6483b4
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/geometry.pro
@@ -0,0 +1,64 @@
+# This file should be viewed as a -*- mode: Makefile -*-
+
+# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
+#			      W A R N I N G ! ! !                             #
+#             a u t h o r i z e d    p e r s o n a l    o n l y               #
+# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
+
+include(../Config.pri)
+
+TEMPLATE        = lib
+TARGET          = $${LIB_GEOMETRY}
+VERSION         = $${APPVERSION}
+TARGET_VERSION_EXT = $${APPVERSION_MAJ}.$${APPVERSION_MID}
+
+#
+# CONFIG
+#
+#######################################
+
+CONFIG          *= dll
+#
+# DEFINES
+#
+#######################################
+
+win32:DEFINES   *= MAKE_LIB_GEOMETRY_DLL
+
+#
+# INCLUDE PATH
+#
+#######################################
+
+#INCLUDEPATH     *= ../system
+
+#
+# BUILD DIRECTORIES
+#
+#######################################
+
+BUILD_DIR       = ../../build
+
+OBJECTS_DIR     = $${BUILD_DIR}/$${REL_OBJECTS_DIR}
+!win32:DESTDIR  = $${BUILD_DIR}/$${REL_DESTDIR}/lib
+win32:DESTDIR   = $${BUILD_DIR}/$${REL_DESTDIR}
+
+#
+# INSTALL
+#
+#######################################
+
+LIB_DIR       = ../../lib
+# install library
+target.path   = $$LIB_DIR
+# "make install" configuration options
+INSTALLS     += target
+
+#
+# SOURCES & HEADERS
+#
+#######################################
+
+!static {
+  include(src.pri)
+}
diff --git a/source/blender/freestyle/intern/geometry/matrix_util.cpp b/source/blender/freestyle/intern/geometry/matrix_util.cpp
new file mode 100755
index 00000000000..2117b06e62f
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/matrix_util.cpp
@@ -0,0 +1,265 @@
+/*
+ *  GXML/Graphite: Geometry and Graphics Programming Library + Utilities
+ *  Copyright (C) 2000 Bruno Levy
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ *
+ *  If you modify this software, you should include a notice giving the
+ *  name of the person performing the modification, the date of modification,
+ *  and the reason for such modification.
+ *
+ *  Contact: Bruno Levy
+ *
+ *     levy@loria.fr
+ *
+ *     ISA Project
+ *     LORIA, INRIA Lorraine, 
+ *     Campus Scientifique, BP 239
+ *     54506 VANDOEUVRE LES NANCY CEDEX 
+ *     FRANCE
+ *
+ *  Note that the GNU General Public License does not permit incorporating
+ *  the Software into proprietary programs. 
+ */
+
+
+#include "matrix_util.h"
+#include <math.h>
+
+
+
+namespace OGF {
+    
+    namespace MatrixUtil {
+
+        static const double EPS = 0.00001 ;
+        static int MAX_ITER = 100 ;
+        
+        void semi_definite_symmetric_eigen(
+            const double *mat, int n, double *eigen_vec, double *eigen_val
+        ) {
+            double *a,*v;
+            double a_norm,a_normEPS,thr,thr_nn;
+            int nb_iter = 0;
+            int jj;
+            int i,j,k,ij,ik,l,m,lm,mq,lq,ll,mm,imv,im,iq,ilv,il,nn;
+            int *index;
+            double a_ij,a_lm,a_ll,a_mm,a_im,a_il;
+            double a_lm_2;
+            double v_ilv,v_imv;
+            double x;
+            double sinx,sinx_2,cosx,cosx_2,sincos;
+            double delta;
+        
+            // Number of entries in mat
+        
+            nn = (n*(n+1))/2;
+        
+            // Step 1: Copy mat to a
+        
+            a = new double[nn];
+        
+            for( ij=0; ij<nn; ij++ ) {
+                a[ij] = mat[ij];
+            }
+        
+            // Ugly Fortran-porting trick: indices for a are between 1 and n
+            a--;
+        
+            // Step 2 : Init diagonalization matrix as the unit matrix
+            v = new double[n*n];
+        
+            ij = 0;
+            for( i=0; i<n; i++ ) {
+                for( j=0; j<n; j++ ) {
+                    if( i==j ) {
+                        v[ij++] = 1.0;
+                    } else {
+                        v[ij++] = 0.0;
+                    }
+                }
+            }
+        
+            // Ugly Fortran-porting trick: indices for v are between 1 and n
+            v--;
+        
+            // Step 3 : compute the weight of the non diagonal terms 
+            ij     = 1  ;
+            a_norm = 0.0;
+            for( i=1; i<=n; i++ ) {
+                for( j=1; j<=i; j++ ) {
+                    if( i!=j ) {
+                        a_ij    = a[ij];
+                        a_norm += a_ij*a_ij;
+                    }
+                    ij++;
+                }
+            }
+        
+            if( a_norm != 0.0 ) {
+        
+                a_normEPS = a_norm*EPS;
+                thr       = a_norm    ;
+    
+                // Step 4 : rotations
+                while( thr > a_normEPS   &&   nb_iter < MAX_ITER ) {
+        
+                    nb_iter++;
+                    thr_nn = thr / nn;
+            
+                    for( l=1  ; l< n; l++ ) {
+                        for( m=l+1; m<=n; m++ ) {
+                    
+                            // compute sinx and cosx 
+                    
+                            lq   = (l*l-l)/2;
+                            mq   = (m*m-m)/2;
+                    
+                            lm     = l+mq;
+                            a_lm   = a[lm];
+                            a_lm_2 = a_lm*a_lm;
+                    
+                            if( a_lm_2 < thr_nn ) {
+                                continue ;
+                            }
+                    
+                            ll   = l+lq;
+                            mm   = m+mq;
+                            a_ll = a[ll];
+                            a_mm = a[mm];
+                    
+                            delta = a_ll - a_mm;
+                    
+                            if( delta == 0.0 ) {
+                                x = - M_PI/4 ; 
+                            } else {
+                                x = - atan( (a_lm+a_lm) / delta ) / 2.0 ;
+                            }
+
+                            sinx    = sin(x)   ;
+                            cosx    = cos(x)   ;
+                            sinx_2  = sinx*sinx;
+                            cosx_2  = cosx*cosx;
+                            sincos  = sinx*cosx;
+                    
+                            // rotate L and M columns 
+       
+                            ilv = n*(l-1);
+                            imv = n*(m-1);
+                    
+                            for( i=1; i<=n;i++ ) {
+                                if( (i!=l) && (i!=m) ) {
+                                    iq = (i*i-i)/2;
+                            
+                                    if( i<m )  { 
+                                        im = i + mq; 
+                                    } else {
+                                        im = m + iq;
+                                    }
+                                    a_im = a[im];
+                            
+                                    if( i<l ) {
+                                        il = i + lq; 
+                                    } else {
+                                        il = l + iq;
+                                    }
+                                    a_il = a[il];
+                            
+                                    a[il] =  a_il*cosx - a_im*sinx;
+                                    a[im] =  a_il*sinx + a_im*cosx;
+                                }
+                        
+                                ilv++;
+                                imv++;
+                        
+                                v_ilv = v[ilv];
+                                v_imv = v[imv];
+                        
+                                v[ilv] = cosx*v_ilv - sinx*v_imv;
+                                v[imv] = sinx*v_ilv + cosx*v_imv;
+                            } 
+                    
+                            x = a_lm*sincos; x+=x;
+                    
+                            a[ll] =  a_ll*cosx_2 + a_mm*sinx_2 - x;
+                            a[mm] =  a_ll*sinx_2 + a_mm*cosx_2 + x;
+                            a[lm] =  0.0;
+                    
+                            thr = fabs( thr - a_lm_2 );
+                        }
+                    }
+                }         
+            }
+        
+            // Step 5: index conversion and copy eigen values 
+        
+            // back from Fortran to C++
+            a++;
+        
+            for( i=0; i<n; i++ ) {
+                k = i + (i*(i+1))/2;
+                eigen_val[i] = a[k];
+            }
+     
+            delete[] a;
+        
+            // Step 6: sort the eigen values and eigen vectors 
+        
+            index = new int[n];
+            for( i=0; i<n; i++ ) {
+                index[i] = i;
+            }
+        
+            for( i=0; i<(n-1); i++ ) {
+                x = eigen_val[i];
+                k = i;
+            
+                for( j=i+1; j<n; j++ ) {
+                    if( x < eigen_val[j] ) {
+                        k = j;
+                        x = eigen_val[j];
+                    }
+                }
+            
+                eigen_val[k] = eigen_val[i];
+                eigen_val[i] = x;
+          
+                jj       = index[k];
+                index[k] = index[i];
+                index[i] = jj;
+            }
+
+
+            // Step 7: save the eigen vectors 
+    
+            v++; // back from Fortran to to C++
+        
+            ij = 0;
+            for( k=0; k<n; k++ ) {
+                ik = index[k]*n;
+                for( i=0; i<n; i++ ) {
+                    eigen_vec[ij++] = v[ik++];
+                }
+            }
+    
+            delete[] v    ;
+            delete[] index;
+            return;
+        }
+    
+//_________________________________________________________
+
+    }
+}
diff --git a/source/blender/freestyle/intern/geometry/matrix_util.h b/source/blender/freestyle/intern/geometry/matrix_util.h
new file mode 100755
index 00000000000..a990413c403
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/matrix_util.h
@@ -0,0 +1,69 @@
+/*
+ *  GXML/Graphite: Geometry and Graphics Programming Library + Utilities
+ *  Copyright (C) 2000 Bruno Levy
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ *
+ *  If you modify this software, you should include a notice giving the
+ *  name of the person performing the modification, the date of modification,
+ *  and the reason for such modification.
+ *
+ *  Contact: Bruno Levy
+ *
+ *     levy@loria.fr
+ *
+ *     ISA Project
+ *     LORIA, INRIA Lorraine, 
+ *     Campus Scientifique, BP 239
+ *     54506 VANDOEUVRE LES NANCY CEDEX 
+ *     FRANCE
+ *
+ *  Note that the GNU General Public License does not permit incorporating
+ *  the Software into proprietary programs. 
+ */
+
+#ifndef __MATRIX_UTIL__
+#define __MATRIX_UTIL__
+
+# include "../system/FreestyleConfig.h"
+
+namespace OGF {
+
+    namespace MatrixUtil {
+
+        /**
+         * computes the eigen values and eigen vectors   
+         * of a semi definite symmetric matrix
+         *
+         * @param  matrix is stored in column symmetric storage, i.e.
+         *     matrix = { m11, m12, m22, m13, m23, m33, m14, m24, m34, m44 ... }
+         *     size = n(n+1)/2
+         *
+         * @param eigen_vectors (return) = { v1, v2, v3, ..., vn } 
+         *   where vk = vk0, vk1, ..., vkn
+         *     size = n^2, must be allocated by caller
+         *
+         * @param eigen_values  (return) are in decreasing order
+         *     size = n,   must be allocated by caller
+         */
+        LIB_GEOMETRY_EXPORT 
+        void semi_definite_symmetric_eigen(
+            const double *mat, int n, double *eigen_vec, double *eigen_val
+        ) ;
+    
+    }
+}
+
+#endif
diff --git a/source/blender/freestyle/intern/geometry/normal_cycle.cpp b/source/blender/freestyle/intern/geometry/normal_cycle.cpp
new file mode 100755
index 00000000000..b456ced8331
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/normal_cycle.cpp
@@ -0,0 +1,103 @@
+
+//
+//  Copyright (C) : Please refer to the COPYRIGHT file distributed 
+//   with this source distribution. 
+//
+//  This program is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU General Public License
+//  as published by the Free Software Foundation; either version 2
+//  of the License, or (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+//
+///////////////////////////////////////////////////////////////////////////////
+#include "normal_cycle.h"
+#include "matrix_util.h"
+
+
+namespace OGF {
+
+//_________________________________________________________
+    
+
+    NormalCycle::NormalCycle() {
+    }
+
+    void NormalCycle::begin() {
+        M_[0] = M_[1] = M_[2] = M_[3] = M_[4] = M_[5] = 0 ;
+    }
+    
+    void NormalCycle::end() {
+        
+        double eigen_vectors[9] ;
+        MatrixUtil::semi_definite_symmetric_eigen(M_, 3, eigen_vectors, eigen_value_) ;
+            
+        axis_[0] = Vec3r(
+            eigen_vectors[0], eigen_vectors[1], eigen_vectors[2]
+        ) ;
+            
+        axis_[1] = Vec3r(
+            eigen_vectors[3], eigen_vectors[4], eigen_vectors[5]
+        ) ;
+        
+        axis_[2] = Vec3r(
+            eigen_vectors[6], eigen_vectors[7], eigen_vectors[8]
+        ) ;
+        
+        // Normalize the eigen vectors
+        
+        for(int i=0; i<3; i++) {
+            axis_[i].normalize() ;
+        }
+
+        // Sort the eigen vectors
+
+        i_[0] = 0 ;
+        i_[1] = 1 ;
+        i_[2] = 2 ;
+
+        double l0 = ::fabs(eigen_value_[0]) ;
+        double l1 = ::fabs(eigen_value_[1]) ;
+        double l2 = ::fabs(eigen_value_[2]) ;
+        
+        if(l1 > l0) {
+            ogf_swap(l0   , l1   ) ;
+            ogf_swap(i_[0], i_[1]) ;
+        }
+        if(l2 > l1) {
+            ogf_swap(l1   , l2   ) ;
+            ogf_swap(i_[1], i_[2]) ;
+        }
+        if(l1 > l0) {
+            ogf_swap(l0   , l1  ) ;
+            ogf_swap(i_[0],i_[1]) ;
+        }
+
+    }
+
+    void NormalCycle::accumulate_dihedral_angle(
+        const Vec3r& edge, double beta, double neigh_area
+    ) {
+        Vec3r e = edge ;
+        e.normalize() ;
+
+        double s = edge.norm() * beta * neigh_area ;
+
+        M_[0] += s * e.x() * e.x() ;
+        M_[1] += s * e.x() * e.y() ;
+        M_[2] += s * e.y() * e.y() ;
+        M_[3] += s * e.x() * e.z() ;
+        M_[4] += s * e.y() * e.z() ;
+        M_[5] += s * e.z() * e.z() ;        
+    }
+    
+//_________________________________________________________
+
+}
diff --git a/source/blender/freestyle/intern/geometry/normal_cycle.h b/source/blender/freestyle/intern/geometry/normal_cycle.h
new file mode 100755
index 00000000000..41fbf7b3fab
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/normal_cycle.h
@@ -0,0 +1,97 @@
+/*
+ *  OGF/Graphite: Geometry and Graphics Programming Library + Utilities
+ *  Copyright (C) 2000 Bruno Levy
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ *
+ *  If you modify this software, you should include a notice giving the
+ *  name of the person performing the modification, the date of modification,
+ *  and the reason for such modification.
+ *
+ *  Contact: Bruno Levy
+ *
+ *     levy@loria.fr
+ *
+ *     ISA Project
+ *     LORIA, INRIA Lorraine, 
+ *     Campus Scientifique, BP 239
+ *     54506 VANDOEUVRE LES NANCY CEDEX 
+ *     FRANCE
+ *
+ *  Note that the GNU General Public License does not permit incorporating
+ *  the Software into proprietary programs. 
+ */
+
+#ifndef __MESH_TOOLS_MATH_NORMAL_CYCLE__
+#define __MESH_TOOLS_MATH_NORMAL_CYCLE__
+
+# include "../system/FreestyleConfig.h"
+# include "Geom.h"
+using namespace Geometry;
+ 
+namespace OGF {
+
+template <class T> inline void ogf_swap(T& x, T& y) {
+        T z = x ;
+        x = y ;
+        y = z ;
+    }
+
+//_________________________________________________________
+
+    /**
+     * NormalCycle evaluates the curvature tensor in function
+     * of a set of dihedral angles and associated vectors.
+     * Reference:
+     *    Restricted Delaunay Triangulation and Normal Cycle,
+     *    D. Cohen-Steiner and J.M. Morvan,
+     *    SOCG 2003
+     */
+    class LIB_GEOMETRY_EXPORT NormalCycle {
+    public:
+        NormalCycle() ;
+        void begin() ;
+        void end() ;
+        /**
+         * Note: the specified edge vector needs to be pre-clipped
+         * by the neighborhood.
+         */
+        void accumulate_dihedral_angle(
+            const Vec3r& edge, real angle, real neigh_area = 1.0
+        ) ;
+        const Vec3r& eigen_vector(int i) const { return axis_[i_[i]] ; }
+        real eigen_value(int i) const   { return eigen_value_[i_[i]] ;  } 
+
+        const Vec3r& N() const    { return eigen_vector(2) ; }
+        const Vec3r& Kmax() const { return eigen_vector(1) ; }
+        const Vec3r& Kmin() const { return eigen_vector(0) ; }
+
+        real n() const    { return eigen_value(2) ; }
+        real kmax() const { return eigen_value(1) ; }
+        real kmin() const { return eigen_value(0) ; }
+
+    private:
+        real center_[3] ;
+        Vec3r axis_[3] ;
+        real eigen_value_[3] ;
+        real M_[6] ;
+        int i_[3] ;
+    } ;
+    
+//_________________________________________________________
+
+}
+
+#endif
diff --git a/source/blender/freestyle/intern/geometry/src.pri b/source/blender/freestyle/intern/geometry/src.pri
new file mode 100755
index 00000000000..a35760fe892
--- /dev/null
+++ b/source/blender/freestyle/intern/geometry/src.pri
@@ -0,0 +1,33 @@
+# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
+#			      W A R N I N G ! ! !                             #
+#             a u t h o r i z e d    p e r s o n a l    o n l y               #
+# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
+
+GEOMETRY_DIR = ../geometry
+
+SOURCES *= $${GEOMETRY_DIR}/GeomCleaner.cpp \
+           $${GEOMETRY_DIR}/GeomUtils.cpp \
+           $${GEOMETRY_DIR}/Grid.cpp \
+           $${GEOMETRY_DIR}/FastGrid.cpp \
+           $${GEOMETRY_DIR}/HashGrid.cpp \
+           $${GEOMETRY_DIR}/FitCurve.cpp \
+           $${GEOMETRY_DIR}/Bezier.cpp \
+           $${GEOMETRY_DIR}/Noise.cpp \
+           $${GEOMETRY_DIR}/matrix_util.cpp \
+           $${GEOMETRY_DIR}/normal_cycle.cpp
+
+HEADERS *= $${GEOMETRY_DIR}/BBox.h \
+           $${GEOMETRY_DIR}/FastGrid.h \
+           $${GEOMETRY_DIR}/Geom.h \
+           $${GEOMETRY_DIR}/GeomCleaner.h \
+           $${GEOMETRY_DIR}/GeomUtils.h \
+           $${GEOMETRY_DIR}/Grid.h \
+           $${GEOMETRY_DIR}/HashGrid.h \
+           $${GEOMETRY_DIR}/Polygon.h \
+           $${GEOMETRY_DIR}/SweepLine.h \
+           $${GEOMETRY_DIR}/FitCurve.h \
+           $${GEOMETRY_DIR}/Bezier.h \
+           $${GEOMETRY_DIR}/Noise.h \
+           $${GEOMETRY_DIR}/VecMat.h \
+           $${GEOMETRY_DIR}/matrix_util.h \
+           $${GEOMETRY_DIR}/normal_cycle.h