diff options
author | Maxime Curioni <maxime.curioni@gmail.com> | 2008-05-08 23:16:40 +0400 |
---|---|---|
committer | Maxime Curioni <maxime.curioni@gmail.com> | 2008-05-08 23:16:40 +0400 |
commit | 64e4a3ec9aed6c8abe095e2cd1fe1552f7cde51c (patch) | |
tree | 6c77358bd447b6c2d215324ef48fc12d1f5ae5ca /source/blender/freestyle/intern/geometry | |
parent | cf2e1e2857cfc5b3c2848c7fc6c9d919ac72fabb (diff) | |
parent | 106974a9d2d5caa5188322507980e3d57d2e3517 (diff) |
soc-2008-mxcurioni: merged changes to revision 14747, cosmetic changes for source/blender/freestyle
Diffstat (limited to 'source/blender/freestyle/intern/geometry')
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 |