Optimized view map calculation by Alexander Beels.

* View map calculation has been intensively optimized for speed by
means of:

1) new spatial grid data structures (SphericalGrid for perspective
cameras and BoxGrid for orthographic cameras; automatically switched
based on the camera type);

2) a heuristic grid density calculation algorithm; and

3) new line visibility computation algorithms: A "traditional"
algorithm for emulating old visibility algorithms, and a "cumulative"
algorithm for improved, more consistent line visibility, both exploiting
the new spatial grid data structures for fast ray casting.

A new option "Raycasting Algorithm" was added to allow users to choose
a ray casting (line visibility) algorithm.  Available choices are:

- Normal Ray Casting
- Fast Ray Casting
- Very Fast Ray Casting
- Culled Traditional Visibility Detection
- Unculled Traditional Visibility Detection
- Culled Cumulative Visibility Detection
- Unculled Cumulative Visibility Detection

The first three algorithms are those available in the original
Freestyle (the "normal" ray casting was used unconditionally, though).
The "fast" and "very fast" ray casting algorithms achieve a faster
calculation at the cost of less visibility accuracy.

The last four are newly introduced optimized options.  The culled
versions of the new algorithms will exclude from visibility
calculation those faces that lay outside the camera, which leads to a
faster view map construction.  The unculled counterparts will take all
faces into account.  The unculled visibility algorithms are useful
when culling affects stroke chaining.

The recommended options for users are the culled/unculled cumulative
visibility algorithms.  These options are meant to replace the old
algorithms in the future.

Performance improvements over the old algorithms depend on the scenes
to be rendered.

* Silhouette detection has also been considerably optimized for speed.

Performance gains by this optimization do not depend on scenes.

* Improper handling of error conditions in the view map construction
was fixed.

23 files changed, 3873 insertions, 60 deletions

diff --git a/source/blender/freestyle/intern/view_map/ArbitraryGridDensityProvider.cpp b/source/blender/freestyle/intern/view_map/ArbitraryGridDensityProvider.cpp
new file mode 100644
index 00000000000..15b2b3343cc
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/ArbitraryGridDensityProvider.cpp
@@ -0,0 +1,108 @@
+//
+//  Filename         : ArbitraryGridDensityProvider.cpp
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2011-2-5
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 "ArbitraryGridDensityProvider.h"
+
+ArbitraryGridDensityProvider::ArbitraryGridDensityProvider(OccluderSource& source, const real proscenium[4], unsigned numCells) 
+	: GridDensityProvider(source), numCells(numCells)
+{
+	initialize (proscenium);
+}
+
+ArbitraryGridDensityProvider::ArbitraryGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform, unsigned numCells) 
+	: GridDensityProvider(source), numCells(numCells)
+{
+	real proscenium[4];
+	calculateQuickProscenium(transform, bbox, proscenium);
+	
+	initialize (proscenium);
+}
+
+ArbitraryGridDensityProvider::ArbitraryGridDensityProvider(OccluderSource& source, unsigned numCells) 
+	: GridDensityProvider(source), numCells(numCells) 
+{
+	real proscenium[4];
+	calculateOptimalProscenium(source, proscenium);
+
+	initialize (proscenium);
+}
+
+ArbitraryGridDensityProvider::~ArbitraryGridDensityProvider () {}
+
+void ArbitraryGridDensityProvider::initialize (const real proscenium[4]) 
+{
+	float prosceniumWidth = (proscenium[1] - proscenium[0]);
+	float prosceniumHeight = (proscenium[3] - proscenium[2]);
+	real cellArea = prosceniumWidth * prosceniumHeight / numCells;
+	cout << prosceniumWidth << " x " << prosceniumHeight << " grid with cells of area " << cellArea << "." << endl;
+
+	_cellSize = sqrt(cellArea);
+	// Now we know how many cells make each side of our grid
+	_cellsX = ceil(prosceniumWidth / _cellSize);
+	_cellsY = ceil(prosceniumHeight / _cellSize);
+	cout << _cellsX << "x" << _cellsY << " cells of size " << _cellSize << " square." << endl;
+
+	// Make sure the grid exceeds the proscenium by a small amount
+	float safetyZone = 0.1;
+	if ( _cellsX * _cellSize < prosceniumWidth * (1.0 + safetyZone) ) {
+		_cellsX = prosceniumWidth * (1.0 + safetyZone) / _cellSize;
+	}
+	if ( _cellsY * _cellSize < prosceniumHeight * (1.0 + safetyZone) ) {
+		_cellsY = prosceniumHeight * (1.0 + safetyZone) / _cellSize;
+	}
+	cout << _cellsX << "x" << _cellsY << " cells of size " << _cellSize << " square." << endl;
+
+	// Find grid origin
+	_cellOrigin[0] = ((proscenium[0] + proscenium[1]) / 2.0) - (_cellsX / 2.0) * _cellSize;
+	_cellOrigin[1] = ((proscenium[2] + proscenium[3]) / 2.0) - (_cellsY / 2.0) * _cellSize;
+}
+
+ArbitraryGridDensityProviderFactory::ArbitraryGridDensityProviderFactory(unsigned numCells)
+	: numCells(numCells)
+{
+}
+
+ArbitraryGridDensityProviderFactory::~ArbitraryGridDensityProviderFactory () {}
+
+auto_ptr<GridDensityProvider> ArbitraryGridDensityProviderFactory::newGridDensityProvider(OccluderSource& source, const real proscenium[4]) 
+{
+	return auto_ptr<GridDensityProvider>(new ArbitraryGridDensityProvider(source, proscenium, numCells));
+}
+
+auto_ptr<GridDensityProvider> ArbitraryGridDensityProviderFactory::newGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform) 
+{
+	return auto_ptr<GridDensityProvider>(new ArbitraryGridDensityProvider(source, bbox, transform, numCells));
+}
+
+auto_ptr<GridDensityProvider> ArbitraryGridDensityProviderFactory::newGridDensityProvider(OccluderSource& source) 
+{
+	return auto_ptr<GridDensityProvider>(new ArbitraryGridDensityProvider(source, numCells));
+}
+
diff --git a/source/blender/freestyle/intern/view_map/ArbitraryGridDensityProvider.h b/source/blender/freestyle/intern/view_map/ArbitraryGridDensityProvider.h
new file mode 100644
index 00000000000..f863b2132a7
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/ArbitraryGridDensityProvider.h
@@ -0,0 +1,67 @@
+//
+//  Filename         : ArbitraryGridDensityProvider.h
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2011-2-5
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 ARBITRARYGRIDDENSITYPROVIDER_H
+#define ARBITRARYGRIDDENSITYPROVIDER_H
+
+#include "GridDensityProvider.h"
+
+class ArbitraryGridDensityProvider : public GridDensityProvider {
+	// Disallow copying and assignment
+	ArbitraryGridDensityProvider (const ArbitraryGridDensityProvider& other);
+	ArbitraryGridDensityProvider& operator= (const ArbitraryGridDensityProvider& other);
+
+public:
+	ArbitraryGridDensityProvider(OccluderSource& source, const real proscenium[4], unsigned numCells);
+	ArbitraryGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform, unsigned numCells);
+	ArbitraryGridDensityProvider(OccluderSource& source, unsigned numCells);
+	virtual ~ArbitraryGridDensityProvider ();
+
+protected:
+	unsigned numCells;
+
+private:
+	void initialize (const real proscenium[4]);
+};
+
+class ArbitraryGridDensityProviderFactory : public GridDensityProviderFactory {
+public:
+	ArbitraryGridDensityProviderFactory(unsigned numCells);
+	~ArbitraryGridDensityProviderFactory ();
+
+	auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source, const real proscenium[4]);
+	auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform);
+	auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source);
+protected:
+	unsigned numCells;
+};
+
+#endif // ARBITRARYGRIDDENSITYPROVIDER_H
+
diff --git a/source/blender/freestyle/intern/view_map/AverageAreaGridDensityProvider.cpp b/source/blender/freestyle/intern/view_map/AverageAreaGridDensityProvider.cpp
new file mode 100644
index 00000000000..8b4c60a7fea
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/AverageAreaGridDensityProvider.cpp
@@ -0,0 +1,120 @@
+//
+//  Filename         : AverageAreaGridDensityProvider.cpp
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2011-2-9
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 "AverageAreaGridDensityProvider.h"
+
+AverageAreaGridDensityProvider::AverageAreaGridDensityProvider(OccluderSource& source, const real proscenium[4], real sizeFactor) 
+	: GridDensityProvider(source)
+{
+	initialize (proscenium, sizeFactor);
+}
+
+AverageAreaGridDensityProvider::AverageAreaGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform, real sizeFactor) 
+	: GridDensityProvider(source)
+{
+	real proscenium[4];
+	calculateQuickProscenium(transform, bbox, proscenium);
+	
+	initialize (proscenium, sizeFactor);
+}
+
+AverageAreaGridDensityProvider::AverageAreaGridDensityProvider(OccluderSource& source, real sizeFactor) 
+	: GridDensityProvider(source)
+{
+	real proscenium[4];
+	calculateOptimalProscenium(source, proscenium);
+
+	initialize (proscenium, sizeFactor);
+}
+
+AverageAreaGridDensityProvider::~AverageAreaGridDensityProvider () {}
+
+void AverageAreaGridDensityProvider::initialize (const real proscenium[4], real sizeFactor) 
+{
+	float prosceniumWidth = (proscenium[1] - proscenium[0]);
+	float prosceniumHeight = (proscenium[3] - proscenium[2]);
+
+	real cellArea = 0.0;
+	unsigned numFaces = 0;
+	for ( source.begin(); source.isValid(); source.next() ) {
+		Polygon3r& poly(source.getGridSpacePolygon());
+		Vec3r min, max;
+		poly.getBBox(min, max);
+		cellArea += (max[0] - min[0]) * (max[1] - min[1]);
+		++numFaces;
+	}
+	cout << "Total area: " << cellArea << ".  Number of faces: " << numFaces << "." << endl;
+	cellArea /= numFaces;
+	cellArea *= sizeFactor;
+	cout << "Building grid with average area " << cellArea << endl;
+
+	_cellSize = sqrt(cellArea);
+	// Now we know how many cells make each side of our grid
+	_cellsX = ceil(prosceniumWidth / _cellSize);
+	_cellsY = ceil(prosceniumHeight / _cellSize);
+	cout << _cellsX << "x" << _cellsY << " cells of size " << _cellSize << " square." << endl;
+
+	// Make sure the grid exceeds the proscenium by a small amount
+	float safetyZone = 0.1;
+	if ( _cellsX * _cellSize < prosceniumWidth * (1.0 + safetyZone) ) {
+		_cellsX = prosceniumWidth * (1.0 + safetyZone) / _cellSize;
+	}
+	if ( _cellsY * _cellSize < prosceniumHeight * (1.0 + safetyZone) ) {
+		_cellsY = prosceniumHeight * (1.0 + safetyZone) / _cellSize;
+	}
+	cout << _cellsX << "x" << _cellsY << " cells of size " << _cellSize << " square." << endl;
+
+	// Find grid origin
+	_cellOrigin[0] = ((proscenium[0] + proscenium[1]) / 2.0) - (_cellsX / 2.0) * _cellSize;
+	_cellOrigin[1] = ((proscenium[2] + proscenium[3]) / 2.0) - (_cellsY / 2.0) * _cellSize;
+}
+
+AverageAreaGridDensityProviderFactory::AverageAreaGridDensityProviderFactory(real sizeFactor)
+	: sizeFactor(sizeFactor)
+{
+}
+
+AverageAreaGridDensityProviderFactory::~AverageAreaGridDensityProviderFactory () {}
+
+auto_ptr<GridDensityProvider> AverageAreaGridDensityProviderFactory::newGridDensityProvider(OccluderSource& source, const real proscenium[4]) 
+{
+	return auto_ptr<GridDensityProvider>(new AverageAreaGridDensityProvider(source, proscenium, sizeFactor));
+}
+
+auto_ptr<GridDensityProvider> AverageAreaGridDensityProviderFactory::newGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform) 
+{
+	return auto_ptr<GridDensityProvider>(new AverageAreaGridDensityProvider(source, bbox, transform, sizeFactor));
+}
+
+auto_ptr<GridDensityProvider> AverageAreaGridDensityProviderFactory::newGridDensityProvider(OccluderSource& source) 
+{
+	return auto_ptr<GridDensityProvider>(new AverageAreaGridDensityProvider(source, sizeFactor));
+}
+
diff --git a/source/blender/freestyle/intern/view_map/AverageAreaGridDensityProvider.h b/source/blender/freestyle/intern/view_map/AverageAreaGridDensityProvider.h
new file mode 100644
index 00000000000..73d28f006a7
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/AverageAreaGridDensityProvider.h
@@ -0,0 +1,67 @@
+//
+//  Filename         : AverageAreaGridDensityProvider.h
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2011-2-9
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 AVERAGEAREAGRIDDENSITYPROVIDER_H
+#define AVERAGEAREAGRIDDENSITYPROVIDER_H
+
+#include "GridDensityProvider.h"
+
+class AverageAreaGridDensityProvider : public GridDensityProvider {
+	// Disallow copying and assignment
+	AverageAreaGridDensityProvider (const AverageAreaGridDensityProvider& other);
+	AverageAreaGridDensityProvider& operator= (const AverageAreaGridDensityProvider& other);
+
+public:
+	AverageAreaGridDensityProvider(OccluderSource& source, const real proscenium[4], real sizeFactor);
+	AverageAreaGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform, real sizeFactor);
+	AverageAreaGridDensityProvider(OccluderSource& source, real sizeFactor);
+	virtual ~AverageAreaGridDensityProvider ();
+
+protected:
+
+private:
+	void initialize (const real proscenium[4], real sizeFactor); 
+};
+
+class AverageAreaGridDensityProviderFactory : public GridDensityProviderFactory {
+public:
+	AverageAreaGridDensityProviderFactory(real sizeFactor);
+	~AverageAreaGridDensityProviderFactory ();
+
+	auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source, const real proscenium[4]);
+	auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform);
+	auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source);
+
+protected:
+	real sizeFactor;
+};
+
+#endif // AVERAGEAREAGRIDDENSITYPROVIDER_H
+
diff --git a/source/blender/freestyle/intern/view_map/BoxGrid.cpp b/source/blender/freestyle/intern/view_map/BoxGrid.cpp
new file mode 100644
index 00000000000..757cf7b6559
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/BoxGrid.cpp
@@ -0,0 +1,210 @@
+//
+//  Filename         : BoxGrid.cpp
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2011-1-29
+//
+///////////////////////////////////////////////////////////////////////////////
+
+//
+//  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 "BoxGrid.h"
+
+#include <stdexcept>
+#include <algorithm>
+
+using namespace std;
+
+// Helper Classes
+
+// OccluderData
+///////////////
+
+// Cell
+/////////
+
+BoxGrid::Cell::Cell () {}
+
+BoxGrid::Cell::~Cell () {}
+
+void BoxGrid::Cell::setDimensions(real x, real y, real sizeX, real sizeY) {
+	const real epsilon = 1.0e-06;
+	boundary[0] = x - epsilon;
+	boundary[1] = x + sizeX + epsilon;
+	boundary[2] = y - epsilon;
+	boundary[3] = y + sizeY + epsilon;
+}
+
+bool BoxGrid::Cell::compareOccludersByShallowestPoint (const BoxGrid::OccluderData* a, const BoxGrid::OccluderData* b) {
+	return a->shallowest < b->shallowest;
+}
+
+void BoxGrid::Cell::indexPolygons() {
+	// Sort occluders by their shallowest points.
+	sort(faces.begin(), faces.end(), compareOccludersByShallowestPoint);
+}
+
+// Iterator
+//////////////////
+
+BoxGrid::Iterator::Iterator (BoxGrid& grid, Vec3r& center, real epsilon) 
+	: _target(grid.transform(center)), 
+	_foundOccludee(false)
+{
+	// Find target cell
+	_cell = grid.findCell(_target);
+	#if boxgridlogging == 1
+		cout << "Searching for occluders of edge centered at " << _target << " in cell [" 
+			<< _cell->boundary[0] << ", " << _cell->boundary[1] << ", " << _cell->boundary[2] 
+			<< ", " << _cell->boundary[3] << "] (" << _cell->faces.size() << " occluders)" << endl;
+	#endif
+
+	// Set iterator
+	_current = _cell->faces.begin();
+}
+
+BoxGrid::Iterator::~Iterator () {}
+
+// BoxGrid
+/////////////////
+
+BoxGrid::BoxGrid(OccluderSource& source, GridDensityProvider& density, ViewMap *viewMap, Vec3r& viewpoint, bool enableQI)
+	: _viewpoint(viewpoint),
+	_enableQI(enableQI)
+{
+	cout << "Generate Cell structure" << endl;
+	// Generate Cell structure
+	assignCells(source, density, viewMap);
+	cout << "Distribute occluders" << endl;
+	// Fill Cells
+	distributePolygons(source);
+	cout << "Reorganize cells" << endl;
+	// Reorganize Cells
+	reorganizeCells();
+	cout << "Ready to use BoxGrid" << endl;
+}
+
+BoxGrid::~BoxGrid () {
+}
+
+void BoxGrid::assignCells (OccluderSource& source, GridDensityProvider& density, ViewMap *viewMap) {
+	_cellSize = density.cellSize();
+	_cellsX = density.cellsX();
+	_cellsY = density.cellsY();
+	_cellOrigin[0] = density.cellOrigin(0);
+	_cellOrigin[1] = density.cellOrigin(1);
+
+	// Now allocate the cell table and fill it with default (empty) cells
+	_cells.resize(_cellsX * _cellsY);
+	for ( cellContainer::iterator i = _cells.begin(), end = _cells.end(); i != end; ++i ) {
+		(*i) = NULL;
+	}
+
+	// Identify cells that will be used, and set the dimensions for each
+	ViewMap::fedges_container& fedges = viewMap->FEdges();
+	for (ViewMap::fedges_container::iterator f = fedges.begin(), fend = fedges.end(); f != fend; ++f ) {
+		if ( (*f)->isInImage() ) {
+			Vec3r point = transform((*f)->center3d());
+			unsigned i, j;
+			getCellCoordinates(point, i, j);
+			if ( _cells[i * _cellsY + j] == NULL ) {
+				// This is an uninitialized cell
+
+				real x, y, width, height;
+
+				x = _cellOrigin[0] + _cellSize * i;
+				width = _cellSize;
+
+				y = _cellOrigin[1] + _cellSize * j;
+				height = _cellSize;
+
+				// Initialize cell
+				Cell* b = _cells[i * _cellsY + j] = new Cell();
+				b->setDimensions(x, y, width, height);
+			}
+		}
+	}
+}
+
+void BoxGrid::distributePolygons (OccluderSource& source) {
+	unsigned long nFaces = 0;
+	unsigned long nKeptFaces = 0;
+
+	for ( source.begin(); source.isValid(); source.next() ) {
+		OccluderData* occluder = NULL;
+
+		try {
+			if ( insertOccluder(source, occluder) ) {
+				_faces.push_back(occluder);
+				++nKeptFaces;
+			}
+		} catch (...) {
+			// If an exception was thrown, _faces.push_back() cannot have succeeded.
+			// occluder is not owned by anyone, and must be deleted.
+			// If the exception was thrown before or during new OccluderData(), then
+			// occluder is NULL, and this delete is harmless.
+			delete occluder;
+			throw;
+		}
+		++nFaces;
+	}
+	cout << "Distributed " << nFaces << " occluders.  Retained " << nKeptFaces << "." << endl;
+}
+
+void BoxGrid::reorganizeCells () {
+	// Sort the occluders by shallowest point
+	for ( vector<Cell*>::iterator i = _cells.begin(), end = _cells.end(); i != end; ++i ) {
+		if ( *i != NULL ) {
+			(*i)->indexPolygons();
+		}
+	}
+}
+
+void BoxGrid::getCellCoordinates(const Vec3r& point, unsigned& x, unsigned& y) {
+	x = min(_cellsX - 1, (unsigned) floor (max((double) 0.0f, point[0] - _cellOrigin[0]) / _cellSize));
+	y = min(_cellsY - 1, (unsigned) floor (max((double) 0.0f, point[1] - _cellOrigin[1]) / _cellSize));
+}
+
+BoxGrid::Cell* BoxGrid::findCell(const Vec3r& point) {
+	unsigned x, y;
+	getCellCoordinates(point, x, y);
+	return _cells[x * _cellsY + y];
+}
+
+bool BoxGrid::orthographicProjection () const {
+	return true;
+}
+
+const Vec3r& BoxGrid::viewpoint() const {
+	return _viewpoint;
+}
+
+bool BoxGrid::enableQI() const {
+	return _enableQI;
+}
+
+BoxGrid::Transform::Transform () : GridHelpers::Transform() {}
+
+Vec3r BoxGrid::Transform::operator() (const Vec3r& point) const {
+	return Vec3r(point[0], point[1], -point[2]);
+}
+
diff --git a/source/blender/freestyle/intern/view_map/BoxGrid.h b/source/blender/freestyle/intern/view_map/BoxGrid.h
new file mode 100644
index 00000000000..43de8d713d5
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/BoxGrid.h
@@ -0,0 +1,376 @@
+//
+//  Filename         : BoxGrid.h
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2011-1-29
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 BOXGRID_H
+#define BOXGRID_H
+
+#define boxgridlogging 0
+
+// I would like to avoid using deque because including ViewMap.h and <deque> or <vector>
+// separately results in redefinitions of identifiers.  ViewMap.h already includes <vector> 
+// so it should be a safe fall-back.
+//#include <vector>
+//#include <deque>
+#include "ViewMap.h"
+#include "../winged_edge/WEdge.h"
+#include "../geometry/Polygon.h"
+#include "../system/PointerSequence.h"
+#include "../geometry/BBox.h"
+#include "../geometry/GridHelpers.h"
+#include "OccluderSource.h"
+#include "GridDensityProvider.h"
+
+class BoxGrid
+{
+public:
+	// Helper classes
+	struct OccluderData {
+		explicit OccluderData (OccluderSource& source, Polygon3r& p);
+		Polygon3r poly;
+		Polygon3r cameraSpacePolygon;
+		real shallowest, deepest;
+		// N.B. We could, of course, store face in poly's userdata
+		// member, like the old ViewMapBuilder code does.  However,
+		// code comments make it clear that userdata is deprecated,
+		// so we avoid the temptation to save 4 or 8 bytes.
+		WFace* face;
+	};
+
+private:
+	struct Cell {
+		// Can't store Cell in a vector without copy and assign
+		//Cell(const Cell& other);
+		//Cell& operator= (const Cell& other);
+
+		explicit Cell ();
+		~Cell ();
+
+		static bool compareOccludersByShallowestPoint (const OccluderData* a, const OccluderData* b);
+
+		void setDimensions(real x, real y, real sizeX, real sizeY);
+		void checkAndInsert(OccluderSource& source, Polygon3r& poly, OccluderData*& occluder);
+		void indexPolygons();
+
+		real boundary[4];
+		//deque<OccluderData*> faces;
+		vector<OccluderData*> faces;
+	};
+
+public:
+	/*****
+
+	Iterator needs to allow the user to avoid full 3D comparison in
+	two cases:
+
+	(1) Where (*current)->deepest < target[2], where the occluder is 
+	unambiguously in front of the target point.
+
+	(2) Where (*current)->shallowest > target[2], where the occluder
+	is unambiguously in back of the target point.
+
+	In addition, when used by OptimizedFindOccludee, Iterator should
+	stop iterating as soon as it has an occludee candidate and 
+	(*current)->shallowest > candidate[2], because at that point forward
+	no new occluder could possibly be a better occludee.
+
+	*****/
+
+	class Iterator {
+	public:
+		// epsilon is not used in this class, but other grids with the same interface may need an epsilon
+		explicit Iterator (BoxGrid& grid, Vec3r& center, real epsilon=1e-06);
+		~Iterator ();
+		void initBeforeTarget ();
+		void initAfterTarget ();
+		void nextOccluder ();
+		void nextOccludee ();
+		bool validBeforeTarget();
+		bool validAfterTarget();
+		WFace* getWFace() const;
+		Polygon3r* getCameraSpacePolygon();
+		void reportDepth(Vec3r origin, Vec3r u, real t);
+	private:
+		bool testOccluder(bool wantOccludee);
+		void markCurrentOccludeeCandidate(real depth);
+
+		Cell* _cell;
+		Vec3r _target;
+		bool _foundOccludee;
+		real _occludeeDepth;
+		//deque<OccluderData*>::iterator _current, _occludeeCandidate;
+		vector<OccluderData*>::iterator _current, _occludeeCandidate;
+	};
+
+	class Transform : public GridHelpers::Transform {
+	public:
+		explicit Transform ();
+		explicit Transform (Transform& other);
+		Vec3r operator() (const Vec3r& point) const;
+	};
+
+private:
+	// Prevent implicit copies and assignments.
+	BoxGrid(const BoxGrid& other);
+	BoxGrid& operator= (const BoxGrid& other);
+public:
+	explicit BoxGrid (OccluderSource& source, GridDensityProvider& density, ViewMap *viewMap, Vec3r& viewpoint, bool enableQI);
+	virtual ~BoxGrid();
+
+	// Generate Cell structure
+	void assignCells(OccluderSource& source, GridDensityProvider& density, ViewMap *viewMap);
+	// Fill Cells
+	void distributePolygons(OccluderSource& source);
+	// Insert one polygon into each matching cell, 
+	// return true if any cell consumes the polygon
+	bool insertOccluder(OccluderSource& source, OccluderData*& occluder);
+	// Sort occluders in each cell
+	void reorganizeCells();
+
+	Cell* findCell(const Vec3r& point);
+
+	// Accessors:
+	bool orthographicProjection() const;
+	const Vec3r& viewpoint() const;
+	bool enableQI() const; 
+	Transform transform;
+
+private:
+	void getCellCoordinates(const Vec3r& point, unsigned& x, unsigned& y);
+
+	typedef PointerSequence<vector<Cell*>, Cell*> cellContainer;
+	//typedef PointerSequence<deque<OccluderData*>, OccluderData*> occluderContainer;
+	typedef PointerSequence<vector<OccluderData*>, OccluderData*> occluderContainer;
+	unsigned _cellsX, _cellsY;
+	float _cellSize;
+	float _cellOrigin[2];
+	cellContainer _cells;
+	occluderContainer _faces;
+	Vec3r _viewpoint;
+	bool _enableQI;
+};
+
+inline void BoxGrid::Iterator::initBeforeTarget () {
+	_current = _cell->faces.begin();
+	while ( _current != _cell->faces.end() && ! testOccluder(false) ) {
+		++_current;
+	}
+}
+
+inline void BoxGrid::Iterator::initAfterTarget () {
+	if ( _foundOccludee ) {
+		#if boxgridlogging == 1
+			std::cout << "\tStarting occludee search from occludeeCandidate at depth " << _occludeeDepth << std::endl;
+		#endif
+		_current = _occludeeCandidate;
+		return;
+	}
+
+	#if boxgridlogging == 1
+		std::cout << "\tStarting occludee search from current position" << std::endl;
+	#endif
+
+	while ( _current != _cell->faces.end() && ! testOccluder(true) ) {
+		++_current;
+	}
+}
+
+inline bool BoxGrid::Iterator::testOccluder (bool wantOccludee) {
+	// End-of-list is not even a valid iterator position
+	if ( _current == _cell->faces.end() ) {
+		// Returning true seems strange, but it will break us out of whatever loop
+		// is calling testOccluder, and _current=_cell->face.end() will make
+		// the calling routine give up.
+		return true;
+	}
+	#if boxgridlogging == 1
+		std::cout << "\tTesting occluder " << (*_current)->poly.getVertices()[0];
+		for ( unsigned i = 1; i < (*_current)->poly.getVertices().size(); ++i ) {
+			std::cout << ", " << (*_current)->poly.getVertices()[i];
+		}
+		std::cout << " from shape " << (*_current)->face->GetVertex(0)->shape()->GetId() << std::endl;
+	#endif
+
+
+	// If we have an occluder candidate and we are unambiguously after it, abort
+	if ( _foundOccludee && (*_current)->shallowest > _occludeeDepth ) {
+		#if boxgridlogging == 1
+			std::cout << "\t\tAborting: shallowest > occludeeCandidate->deepest" << std::endl;
+		#endif
+		_current = _cell->faces.end();
+
+		// See note above
+		return true;
+	}
+
+	// Specific continue or stop conditions when searching for each type
+	if ( wantOccludee ) {
+		if ( (*_current)->deepest < _target[2] ) {
+			#if boxgridlogging == 1
+				std::cout << "\t\tSkipping: shallower than target while looking for occludee" << std::endl;
+			#endif
+			return false;
+		}
+	} else {
+		if ( (*_current)->shallowest > _target[2] ) {
+			#if boxgridlogging == 1
+				std::cout << "\t\tStopping: deeper than target while looking for occluder" << std::endl;
+			#endif
+			return true;
+		}
+	}
+
+	// Depthwise, this is a valid occluder.
+
+	// Check to see if target is in the 2D bounding box
+	Vec3r bbMin, bbMax;
+	(*_current)->poly.getBBox(bbMin, bbMax);
+	if ( _target[0] < bbMin[0] || _target[0] > bbMax[0] || _target[1] < bbMin[1] || _target[1] > bbMax[1] ) {
+		#if boxgridlogging == 1
+			std::cout << "\t\tSkipping: bounding box violation" << std::endl;
+		#endif
+		return false;
+	}
+
+	// We've done all the corner cutting we can.
+	// Let the caller work out whether or not
+	// the geometry is correct.
+	return true;
+}
+
+inline void BoxGrid::Iterator::reportDepth (Vec3r origin, Vec3r u, real t) {
+	// The reported depth is the length of a ray in camera space
+	// We need to convert it into a Z-value in grid space
+	real depth = -(origin + (u * t))[2];
+	#if boxgridlogging == 1
+		std::cout << "\t\tReporting depth of occluder/ee: " << depth;
+	#endif
+	if ( depth > _target[2] ) {
+		#if boxgridlogging == 1
+			std::cout << " is deeper than target" << std::endl;
+		#endif
+		// If the current occluder is the best occludee so far, save it.
+		if ( ! _foundOccludee || _occludeeDepth > depth ) {
+			markCurrentOccludeeCandidate(depth);
+		} 
+	} else {
+		#if boxgridlogging == 1
+			std::cout << std::endl;
+		#endif
+	}
+}
+
+inline void BoxGrid::Iterator::nextOccluder () {
+	if ( _current != _cell->faces.end() ) {
+		do {
+			++_current;
+		} while ( _current != _cell->faces.end() && ! testOccluder(false) );
+	}
+}
+
+inline void BoxGrid::Iterator::nextOccludee () {
+	if ( _current != _cell->faces.end() ) {
+		do {
+			++_current;
+		} while ( _current != _cell->faces.end() && ! testOccluder(true) );
+	}
+}
+
+inline bool BoxGrid::Iterator::validBeforeTarget () {
+	return _current != _cell->faces.end() && (*_current)->shallowest <= _target[2];
+}
+
+inline bool BoxGrid::Iterator::validAfterTarget () {
+	return _current != _cell->faces.end();
+}
+
+inline void BoxGrid::Iterator::markCurrentOccludeeCandidate(real depth) {
+	#if boxgridlogging == 1
+		std::cout << "\t\tFound occludeeCandidate at depth " << depth << std::endl;
+	#endif
+	_occludeeCandidate = _current;
+	_occludeeDepth = depth;
+	_foundOccludee = true;
+}
+
+inline WFace* BoxGrid::Iterator::getWFace() const {
+	return (*_current)->face;
+}
+
+inline Polygon3r* BoxGrid::Iterator::getCameraSpacePolygon() {
+	return &((*_current)->cameraSpacePolygon);
+}
+
+inline BoxGrid::OccluderData::OccluderData (OccluderSource& source, Polygon3r& p)
+	: poly(p), 
+	cameraSpacePolygon(source.getCameraSpacePolygon()),
+	face(source.getWFace())
+{
+	// Set shallowest and deepest based on bbox
+	Vec3r min, max;
+	poly.getBBox(min, max);
+	shallowest = min[2];
+	deepest = max[2];
+}
+
+inline void BoxGrid::Cell::checkAndInsert(OccluderSource& source, Polygon3r& poly, OccluderData*& occluder) {
+	if ( GridHelpers::insideProscenium (boundary, poly) ) {
+		if ( occluder == NULL) {
+			// Disposal of occluder will be handled in BoxGrid::distributePolygons(),
+			// or automatically by BoxGrid::_faces;
+			occluder = new OccluderData(source, poly);
+		}
+		faces.push_back(occluder);
+	}
+}
+
+inline bool BoxGrid::insertOccluder(OccluderSource& source, OccluderData*& occluder) {
+	Polygon3r& poly(source.getGridSpacePolygon());
+	occluder = NULL;
+
+	Vec3r bbMin, bbMax;
+	poly.getBBox(bbMin, bbMax);
+	// Check overlapping cells
+	unsigned startX, startY, endX, endY;
+	getCellCoordinates(bbMin, startX, startY);
+	getCellCoordinates(bbMax, endX, endY);
+
+	for ( unsigned i = startX; i <= endX; ++i ) {
+		for ( unsigned j = startY; j <= endY; ++j ) {
+			if ( _cells[i * _cellsY + j] != NULL ) {
+				_cells[i * _cellsY + j]->checkAndInsert(source, poly, occluder);
+			}
+		}
+	}
+
+	return occluder != NULL;
+}
+
+#endif // BOXGRID_H
+
diff --git a/source/blender/freestyle/intern/view_map/CulledOccluderSource.cpp b/source/blender/freestyle/intern/view_map/CulledOccluderSource.cpp
new file mode 100644
index 00000000000..ea57da93347
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/CulledOccluderSource.cpp
@@ -0,0 +1,277 @@
+//
+//  Filename         : CulledOccluderSource.h
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2010-12-21
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 "CulledOccluderSource.h"
+#include "../geometry/GridHelpers.h"
+#include "FRS_freestyle.h"
+
+CulledOccluderSource::CulledOccluderSource (const GridHelpers::Transform& t, WingedEdge& we, ViewMap& viewMap, bool extensiveFEdgeSearch) 
+	: OccluderSource(t, we), 
+	rejected(0),
+	gridSpaceOccluderProsceniumInitialized(false)
+{
+	cullViewEdges(viewMap, extensiveFEdgeSearch);
+
+	// If we have not found any visible FEdges during our cull, then there is nothing
+	// to iterate over.  Short-circuit everything.
+	valid = gridSpaceOccluderProsceniumInitialized;
+
+	if ( valid && ! testCurrent() ) {
+		next();
+	}
+}
+
+CulledOccluderSource::~CulledOccluderSource() {
+}
+
+bool CulledOccluderSource::testCurrent() {
+	if ( valid ) {
+		// The test for gridSpaceOccluderProsceniumInitialized should not be necessary
+		return gridSpaceOccluderProsceniumInitialized && GridHelpers::insideProscenium (gridSpaceOccluderProscenium, cachedPolygon);
+	}
+	return false;
+}
+
+bool CulledOccluderSource::next() {
+	while ( OccluderSource::next() ) {
+		if ( testCurrent() ) {
+			++rejected;
+			return true;
+		}
+	}
+	std::cout << "Finished generating occluders.  Rejected " << rejected << " faces." << std::endl;
+	return false;
+}
+
+void CulledOccluderSource::getOccluderProscenium(real proscenium[4]) {
+	for ( unsigned i = 0; i < 4; ++i ) {
+		proscenium[i] = gridSpaceOccluderProscenium[i];
+	}
+}
+
+static inline real distance2D(const Vec3r & point, const real origin[2]) {
+	return ::hypot((point[0] - origin[0]), (point[1] - origin[1]));
+}
+
+static inline bool crossesProscenium(real proscenium[4], FEdge *fe) {
+	Vec2r min(proscenium[0], proscenium[2]);
+	Vec2r max(proscenium[1], proscenium[3]);
+	Vec2r A(fe->vertexA()->getProjectedX(), fe->vertexA()->getProjectedY());
+	Vec2r B(fe->vertexB()->getProjectedX(), fe->vertexB()->getProjectedY());
+
+	return GeomUtils::intersect2dSeg2dArea (min, max, A, B);
+}
+
+static inline bool insideProscenium(real proscenium[4], const Vec3r& point) {
+	return ! ( point[0] < proscenium[0] || point[0] > proscenium[1] || point[1] < proscenium[2] || point[1] > proscenium[3] );
+}
+
+void CulledOccluderSource::cullViewEdges(ViewMap& viewMap, bool extensiveFEdgeSearch) {
+	// Cull view edges by marking them as non-displayable.
+	// This avoids the complications of trying to delete
+	// edges from the ViewMap.
+
+	// Non-displayable view edges will be skipped over during
+	// visibility calculation.
+
+	// View edges will be culled according to their position
+	// w.r.t. the viewport proscenium (viewport + 5% border,
+	// or some such).
+
+	// Get proscenium boundary for culling
+	real viewProscenium[4];
+	GridHelpers::getDefaultViewProscenium(viewProscenium);
+	real prosceniumOrigin[2];
+	prosceniumOrigin[0] = (viewProscenium[1] - viewProscenium[0]) / 2.0;
+	prosceniumOrigin[1] = (viewProscenium[3] - viewProscenium[2]) / 2.0;
+	cout << "Proscenium culling:" << endl;
+	cout << "Proscenium: [" << viewProscenium[0] << ", " << viewProscenium[1] << ", " << viewProscenium[2] << ", " << viewProscenium[3] << "]"<< endl;
+	cout << "Origin: [" << prosceniumOrigin[0] << ", " << prosceniumOrigin[1] << "]"<< endl;
+	
+	// A separate occluder proscenium will also be maintained,
+	// starting out the same as the viewport proscenium, and 
+	// expanding as necessary so that it encompasses the center 
+	// point of at least one feature edge in each retained view
+	// edge.
+	// The occluder proscenium will be used later to cull occluding
+	// triangles before they are inserted into the Grid.
+	// The occluder proscenium starts out the same size as the view
+	// proscenium
+	GridHelpers::getDefaultViewProscenium(occluderProscenium);
+
+	// N.B. Freestyle is inconsistent in its use of ViewMap::viewedges_container
+	// and vector<ViewEdge*>::iterator.  Probably all occurences of vector<ViewEdge*>::iterator
+	// should be replaced ViewMap::viewedges_container throughout the code.
+	// For each view edge
+	ViewMap::viewedges_container::iterator ve, veend;
+
+	for(ve=viewMap.ViewEdges().begin(), veend=viewMap.ViewEdges().end(); ve!=veend; ve++) {
+		// Overview:
+		//    Search for a visible feature edge
+		//    If none: mark view edge as non-displayable
+		//    Otherwise:
+		//        Find a feature edge with center point inside occluder proscenium.
+		//        If none exists, find the feature edge with center point
+		//            closest to viewport origin.
+		//            Expand occluder proscenium to enclose center point.
+
+		// For each feature edge, while bestOccluderTarget not found and view edge not visibile
+		bool bestOccluderTargetFound = false;
+		FEdge *bestOccluderTarget = NULL;
+		real bestOccluderDistance = 0.0;
+		FEdge *festart = (*ve)->fedgeA();
+		FEdge *fe = festart;
+		// All ViewEdges start culled
+		(*ve)->setIsInImage(false);
+
+		// For simple visibility calculation: mark a feature edge
+		// that is known to have a center point inside the occluder proscenium.  
+		// Cull all other feature edges.
+		do {
+			// All FEdges start culled
+			fe->setIsInImage(false);
+
+			// Look for the visible edge that can most easily be included
+			// in the occluder proscenium.
+			if ( ! bestOccluderTargetFound ) {
+				// If center point is inside occluder proscenium, 
+				if ( insideProscenium(occluderProscenium, fe->center2d()) ) {
+					// Use this feature edge for visibility deterimination
+					fe->setIsInImage(true);
+					expandGridSpaceOccluderProscenium(fe);
+					// Mark bestOccluderTarget as found
+					bestOccluderTargetFound = true;
+					bestOccluderTarget = fe;
+				} else {
+					real d = distance2D(fe->center2d(), prosceniumOrigin);
+					// If center point is closer to viewport origin than current target
+					if ( bestOccluderTarget == NULL || d < bestOccluderDistance ) {
+						// Then store as bestOccluderTarget
+						bestOccluderDistance = d;
+						bestOccluderTarget = fe;
+					}
+				}
+			}
+
+			// If feature edge crosses the view proscenium
+			if ( ! (*ve)->isInImage() && crossesProscenium(viewProscenium, fe) ) {
+				// Then the view edge will be included in the image
+				(*ve)->setIsInImage(true);
+			}
+			fe = fe->nextEdge();
+		} while ( fe != NULL && fe != festart && ! ( bestOccluderTargetFound && (*ve)->isInImage() ) );
+
+		// Either we have run out of FEdges, or we already have the one edge we need to determine visibility
+		// Cull all remaining edges.
+		while ( fe != NULL && fe != festart ) {
+			fe->setIsInImage(false);
+			fe = fe->nextEdge();
+		}
+
+		// If bestOccluderTarget was not found inside the occluder proscenium, 
+		// we need to expand the occluder proscenium to include it.
+		if ( (*ve)->isInImage() && bestOccluderTarget != NULL && ! bestOccluderTargetFound ) {
+			// Expand occluder proscenium to enclose bestOccluderTarget
+			Vec3r point = bestOccluderTarget->center2d();
+			if ( point[0] < occluderProscenium[0] ) {
+				occluderProscenium[0] = point[0];
+			} else if ( point[0] > occluderProscenium[1] ) {
+				occluderProscenium[1] = point[0];
+			}
+			if ( point[1] < occluderProscenium[2] ) {
+				occluderProscenium[2] = point[1];
+			} else if ( point[1] > occluderProscenium[3] ) {
+				occluderProscenium[3] = point[1];
+			}
+			// Use bestOccluderTarget for visibility determination
+			bestOccluderTarget->setIsInImage(true);
+		}
+	}
+
+	// We are done calculating the occluder proscenium.
+	// Expand the occluder proscenium by an epsilon to avoid rounding errors.
+	const real epsilon = 1.0e-6;
+	occluderProscenium[0] -= epsilon;
+	occluderProscenium[1] += epsilon;
+	occluderProscenium[2] -= epsilon;
+	occluderProscenium[3] += epsilon;
+
+	// For "Normal" or "Fast" style visibility computation only:
+
+	// For more detailed visibility calculation, make a second pass through
+	// the view map, marking all feature edges with center points inside
+	// the final occluder proscenium.  All of these feature edges can be
+	// considered during visibility calculation.
+
+	// So far we have only found one FEdge per ViewEdge.  The "Normal" and
+	// "Fast" styles of visibility computation want to consider many
+	// FEdges for each ViewEdge.
+	// Here we re-scan the view map to find any usable FEdges that we
+	// skipped on the first pass, or that have become usable because the
+	// occluder proscenium has been expanded since the edge was visited
+	// on the first pass.
+	if ( extensiveFEdgeSearch ) {
+		// For each view edge,
+		for(ve=viewMap.ViewEdges().begin(), veend=viewMap.ViewEdges().end(); ve!=veend; ve++) {
+			if ( ! (*ve)->isInImage() ) {
+				continue;
+			}
+			// For each feature edge,
+			FEdge *festart = (*ve)->fedgeA();
+			FEdge *fe = festart;
+			do {
+				// If not (already) visible and center point inside occluder proscenium, 
+				if ( ! fe->isInImage() && insideProscenium(occluderProscenium, fe->center2d()) ) {
+					// Use the feature edge for visibility determination
+					fe->setIsInImage(true);
+					expandGridSpaceOccluderProscenium(fe);
+				}
+				fe = fe->nextEdge();
+			} while ( fe != NULL && fe != festart );
+		}
+	}
+
+	// Up until now, all calculations have been done in camera space.
+	// However, the occluder source's iteration and the grid that consumes the occluders
+	// both work in gridspace, so we need a version of the occluder proscenium in gridspace.
+	// Set the gridspace occlude proscenium
+}
+
+void CulledOccluderSource::expandGridSpaceOccluderProscenium(FEdge* fe) {
+	if ( gridSpaceOccluderProsceniumInitialized ) {
+		GridHelpers::expandProscenium (gridSpaceOccluderProscenium, transform(fe->center3d()));
+	} else {
+		const Vec3r& point = transform(fe->center3d());
+		gridSpaceOccluderProscenium[0] = gridSpaceOccluderProscenium[1] = point[0];
+		gridSpaceOccluderProscenium[2] = gridSpaceOccluderProscenium[3] = point[1];
+		gridSpaceOccluderProsceniumInitialized = true;
+	}
+}
+
diff --git a/source/blender/freestyle/intern/view_map/CulledOccluderSource.h b/source/blender/freestyle/intern/view_map/CulledOccluderSource.h
new file mode 100644
index 00000000000..3c00d5e34ad
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/CulledOccluderSource.h
@@ -0,0 +1,63 @@
+//
+//  Filename         : CulledOccluderSource.h
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2010-12-21
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 CULLEDOCCLUDERSOURCE_H
+#define CULLEDOCCLUDERSOURCE_H
+
+#include "OccluderSource.h"
+#include "ViewMap.h"
+
+class CulledOccluderSource : public OccluderSource {
+	// Disallow copying and assignment
+	CulledOccluderSource (const CulledOccluderSource& other);
+	CulledOccluderSource& operator= (const CulledOccluderSource& other);
+
+public:
+	CulledOccluderSource (const GridHelpers::Transform& transform, WingedEdge& we, ViewMap& viewMap, bool extensiveFEdgeSearch = true);
+	virtual ~CulledOccluderSource();
+
+	void cullViewEdges(ViewMap& viewMap, bool extensiveFEdgeSearch);
+
+	bool next();
+
+	void getOccluderProscenium(real proscenium[4]);
+
+private:
+	bool testCurrent();
+	void expandGridSpaceOccluderProscenium(FEdge* fe);
+
+	real occluderProscenium[4];
+	real gridSpaceOccluderProscenium[4];
+
+	unsigned long rejected;
+	bool gridSpaceOccluderProsceniumInitialized;
+};
+
+#endif // CULLEDOCCLUDERSOURCE_H
diff --git a/source/blender/freestyle/intern/view_map/FEdgeXDetector.cpp b/source/blender/freestyle/intern/view_map/FEdgeXDetector.cpp
index 16c38c63813..65fe146aec5 100755
--- a/source/blender/freestyle/intern/view_map/FEdgeXDetector.cpp
+++ b/source/blender/freestyle/intern/view_map/FEdgeXDetector.cpp
@@ -114,16 +114,18 @@ void FEdgeXDetector::preProcessShape(WXShape* iWShape) {
     preProcessFace((WXFace*)(*f));
   }
 
-  vector<WVertex*>& wvertices = iWShape->getVertexList();
-  for(vector<WVertex*>::iterator wv=wvertices.begin(), wvend=wvertices.end(); 
-  wv!=wvend;
-  ++wv){
-    // Compute curvatures
-    WXVertex * wxv = dynamic_cast<WXVertex*>(*wv);
-    computeCurvatures(wxv);
+  if(_faceSmoothness || _computeRidgesAndValleys || _computeSuggestiveContours ) {
+    vector<WVertex*>& wvertices = iWShape->getVertexList();
+    for(vector<WVertex*>::iterator wv=wvertices.begin(), wvend=wvertices.end(); 
+    wv!=wvend;
+    ++wv){
+      // Compute curvatures
+      WXVertex * wxv = dynamic_cast<WXVertex*>(*wv);
+      computeCurvatures(wxv);
+    }
+    _meanK1 /= (real)(_nPoints);
+    _meanKr /= (real)(_nPoints);
   }
-  _meanK1 /= (real)(_nPoints);
-  _meanKr /= (real)(_nPoints);
 }
 
 void FEdgeXDetector::preProcessFace(WXFace *iFace){
@@ -603,7 +605,7 @@ void FEdgeXDetector::postProcessSuggestiveContourFace(WXFace *iFace) {
   WXVertex *v, *opposite_vertex_a, *opposite_vertex_b;
   WXFace *wxf;
   WOEdge *opposite_edge;
-  Vec3r opposite_edge_vec, normal_vec, radial_normal_vec, er_vec, v_vec, inter, inter1, inter2, tmp_vec;
+  Vec3r normal_vec, radial_normal_vec, er_vec, v_vec, inter, inter1, inter2, tmp_vec;
   GeomUtils::intersection_test res;
   real kr(0), kr1(0), kr2(0), t;
 
@@ -629,12 +631,11 @@ void FEdgeXDetector::postProcessSuggestiveContourFace(WXFace *iFace) {
       
       opposite_vertex_a = (WXVertex*)opposite_edge->GetaVertex();
       opposite_vertex_b = (WXVertex*)opposite_edge->GetbVertex();
-      opposite_edge_vec = opposite_vertex_b->GetVertex() - opposite_vertex_a->GetVertex();
       normal_vec = wxf->GetVertexNormal(v); // FIXME: what about e1 ^ e2 ?
       radial_normal_vec = er_vec ^ normal_vec;
 
       // Test wether the radial plan intersects with the edge at the opposite of v.
-      res = GeomUtils::intersectRayPlane(opposite_vertex_a->GetVertex(), opposite_edge_vec,
+      res = GeomUtils::intersectRayPlane(opposite_vertex_a->GetVertex(), opposite_edge->GetVec(),
 					 radial_normal_vec, -(v_vec * radial_normal_vec),
 					 t,
 					 1.e-06);
@@ -642,7 +643,7 @@ void FEdgeXDetector::postProcessSuggestiveContourFace(WXFace *iFace) {
       // If there is an intersection, compute the value of the derivative ath that point.
       if ((res == GeomUtils::DO_INTERSECT) && (t >= 0) && (t <= 1)) {
 	kr = t * opposite_vertex_a->curvatures()->Kr + (1 - t) * opposite_vertex_b->curvatures()->Kr;
-	inter = opposite_vertex_a->GetVertex() + t * opposite_edge_vec;
+	inter = opposite_vertex_a->GetVertex() + t * opposite_edge->GetVec();
 	tmp_vec = inter - v->GetVertex();
 	// Is it kr1 or kr2?
 	if (tmp_vec * er_vec > 0) {
diff --git a/source/blender/freestyle/intern/view_map/FEdgeXDetector.h b/source/blender/freestyle/intern/view_map/FEdgeXDetector.h
index 4cbaa4a7b5a..d08d75a84b6 100755
--- a/source/blender/freestyle/intern/view_map/FEdgeXDetector.h
+++ b/source/blender/freestyle/intern/view_map/FEdgeXDetector.h
@@ -58,6 +58,7 @@ public:
 	_computeMaterialBoundaries = true;
     _sphereRadius = 1.0;
     _orthographicProjection = false;
+    _faceSmoothness = false;
     _changes = false;
     _kr_derivative_epsilon = 0.0;
 	_creaseAngle = 0.7; // angle of 134.43 degrees
@@ -135,6 +136,12 @@ public:
   inline void enableRidgesAndValleysFlag(bool b) {_computeRidgesAndValleys = b;}
   inline void enableSuggestiveContours(bool b) {_computeSuggestiveContours = b;}
   inline void enableMaterialBoundaries(bool b) {_computeMaterialBoundaries = b;}
+  inline void enableFaceSmoothness(bool b) {
+    if (b != _faceSmoothness) {
+	  _faceSmoothness = b;
+      _changes=true;
+    }
+  }
   /*! Sets the radius of the geodesic sphere around each vertex (for the curvature computation)
    *  \param r
    *    The radius of the sphere expressed as a ratio of the mean edge size
@@ -167,6 +174,7 @@ protected:
   bool _computeRidgesAndValleys;
   bool _computeSuggestiveContours;
   bool _computeMaterialBoundaries;
+  bool _faceSmoothness;
   real _sphereRadius; // expressed as a ratio of the mean edge size
   real _creaseAngle; // [-1, 1] compared with the inner product of face normals
   bool _changes;
diff --git a/source/blender/freestyle/intern/view_map/GridDensityProvider.h b/source/blender/freestyle/intern/view_map/GridDensityProvider.h
new file mode 100644
index 00000000000..078fc5f2c98
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/GridDensityProvider.h
@@ -0,0 +1,131 @@
+//
+//  Filename         : GridDensityProvider.h
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2011-2-5
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 GRIDDENSITYPROVIDER_H
+#define GRIDDENSITYPROVIDER_H
+
+#include <stdexcept>
+#include <memory>
+#include "OccluderSource.h"
+#include "../geometry/BBox.h"
+
+class GridDensityProvider {
+	// Disallow copying and assignment
+	GridDensityProvider (const GridDensityProvider& other);
+	GridDensityProvider& operator= (const GridDensityProvider& other);
+
+public:
+	GridDensityProvider (OccluderSource& source)
+		: source(source) {
+	}
+
+	virtual ~GridDensityProvider() {};
+
+	float cellSize() {
+		return _cellSize;
+	}
+
+	unsigned cellsX() {
+		return _cellsX;
+	}
+
+	unsigned cellsY() {
+		return _cellsY;
+	}
+
+	float cellOrigin(int index) {
+		if ( index < 2 ) {
+			return _cellOrigin[index];
+		} else {
+			throw new out_of_range("GridDensityProvider::cellOrigin can take only indexes of 0 or 1.");
+		}
+	}
+
+	static void calculateOptimalProscenium(OccluderSource& source, real proscenium[4]) {
+		source.begin();
+		if ( source.isValid() ) {
+			const Vec3r& initialPoint = source.getGridSpacePolygon().getVertices()[0];
+			proscenium[0] = proscenium[1] = initialPoint[0];
+			proscenium[2] = proscenium[3] = initialPoint[1];
+			while ( source.isValid() ) {
+				GridHelpers::expandProscenium (proscenium, source.getGridSpacePolygon());
+				source.next();
+			}
+		}
+		cout << "Proscenium: (" << proscenium[0] << ", " << proscenium[1] << ", " << proscenium[2] << ", " << proscenium[3] << ")" << endl;
+	}
+
+	static void calculateQuickProscenium(const GridHelpers::Transform& transform, const BBox<Vec3r>& bbox, real proscenium[4]) {
+		real z;
+		// We want to use the z-coordinate closest to the camera to determine the proscenium face
+		if ( ::fabs(bbox.getMin()[2]) < ::fabs(bbox.getMax()[2]) ) {
+			z = bbox.getMin()[2];
+		} else {
+			z = bbox.getMax()[2];
+		}
+		// Now calculate the proscenium according to the min and max values of the x and y coordinates
+		Vec3r minPoint = transform(Vec3r(bbox.getMin()[0], bbox.getMin()[1], z));
+		Vec3r maxPoint = transform(Vec3r(bbox.getMax()[0], bbox.getMax()[1], z));
+		cout << "Bounding box: " << minPoint << " to " << maxPoint << endl;
+		proscenium[0] = std::min(minPoint[0], maxPoint[0]);
+		proscenium[1] = std::max(minPoint[0], maxPoint[0]);
+		proscenium[2] = std::min(minPoint[1], maxPoint[1]);
+		proscenium[3] = std::max(minPoint[1], maxPoint[1]);
+		cout << "Proscenium  : " << proscenium[0] << ", " << proscenium[1] << ", " << proscenium[2] << ", " << proscenium[3] << endl;
+	}
+
+protected:
+	OccluderSource& source;
+	unsigned _cellsX, _cellsY;
+	float _cellSize;
+	float _cellOrigin[2];
+};
+
+class GridDensityProviderFactory {
+	// Disallow copying and assignment
+	GridDensityProviderFactory (const GridDensityProviderFactory& other);
+	GridDensityProviderFactory& operator= (const GridDensityProviderFactory& other);
+
+public:
+	GridDensityProviderFactory() 
+	{
+	}
+
+	virtual auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source, const real proscenium[4]) =0;
+
+	virtual auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform) =0;
+
+	virtual auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source) =0;
+
+	virtual ~GridDensityProviderFactory () {}
+};
+
+#endif // GRIDDENSITYPROVIDER_H
+
diff --git a/source/blender/freestyle/intern/view_map/HeuristicGridDensityProviderFactory.cpp b/source/blender/freestyle/intern/view_map/HeuristicGridDensityProviderFactory.cpp
new file mode 100644
index 00000000000..05b70af90b3
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/HeuristicGridDensityProviderFactory.cpp
@@ -0,0 +1,74 @@
+//
+//  Filename         : HeuristicGridDensityProviderFactory.cpp
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2011-2-8
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 "HeuristicGridDensityProviderFactory.h"
+
+HeuristicGridDensityProviderFactory::HeuristicGridDensityProviderFactory(real sizeFactor, unsigned numFaces)
+	: sizeFactor(sizeFactor), numFaces(numFaces)
+{
+}
+
+HeuristicGridDensityProviderFactory::~HeuristicGridDensityProviderFactory () {}
+
+auto_ptr<GridDensityProvider> HeuristicGridDensityProviderFactory::newGridDensityProvider(OccluderSource& source, const real proscenium[4]) 
+{
+	auto_ptr<AverageAreaGridDensityProvider> avg(new AverageAreaGridDensityProvider(source, proscenium, sizeFactor));
+	auto_ptr<Pow23GridDensityProvider> p23(new Pow23GridDensityProvider(source, proscenium, numFaces));
+	if ( avg->cellSize() > p23->cellSize() ) {
+		return (auto_ptr<GridDensityProvider>) p23;
+	} else {
+		return (auto_ptr<GridDensityProvider>) avg;
+	}
+}
+
+auto_ptr<GridDensityProvider> HeuristicGridDensityProviderFactory::newGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform) 
+{
+	auto_ptr<AverageAreaGridDensityProvider> avg(new AverageAreaGridDensityProvider(source, bbox, transform, sizeFactor));
+	auto_ptr<Pow23GridDensityProvider> p23(new Pow23GridDensityProvider(source, bbox, transform, numFaces));
+	if ( avg->cellSize() > p23->cellSize() ) {
+		return (auto_ptr<GridDensityProvider>) p23;
+	} else {
+		return (auto_ptr<GridDensityProvider>) avg;
+	}
+}
+
+auto_ptr<GridDensityProvider> HeuristicGridDensityProviderFactory::newGridDensityProvider(OccluderSource& source) 
+{
+	real proscenium[4];
+	GridDensityProvider::calculateOptimalProscenium(source, proscenium);
+	auto_ptr<AverageAreaGridDensityProvider> avg(new AverageAreaGridDensityProvider(source, proscenium, sizeFactor));
+	auto_ptr<Pow23GridDensityProvider> p23(new Pow23GridDensityProvider(source, proscenium, numFaces));
+	if ( avg->cellSize() > p23->cellSize() ) {
+		return (auto_ptr<GridDensityProvider>) p23;
+	} else {
+		return (auto_ptr<GridDensityProvider>) avg;
+	}
+}
+
diff --git a/source/blender/freestyle/intern/view_map/HeuristicGridDensityProviderFactory.h b/source/blender/freestyle/intern/view_map/HeuristicGridDensityProviderFactory.h
new file mode 100644
index 00000000000..d32c4ae4407
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/HeuristicGridDensityProviderFactory.h
@@ -0,0 +1,54 @@
+//
+//  Filename         : HeuristicGridDensityProviderFactory.h
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2011-2-8
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 HEURISTICGRIDDENSITYPROVIDERFACTORY_H
+#define HEURISTICGRIDDENSITYPROVIDERFACTORY_H
+
+// #include <memory> // provided by GridDensityProvider.h
+// #include "GridDensityProvider.h" // provided by *GridDensityProvider.h below
+#include "Pow23GridDensityProvider.h"
+#include "AverageAreaGridDensityProvider.h"
+
+class HeuristicGridDensityProviderFactory : public GridDensityProviderFactory {
+public:
+	HeuristicGridDensityProviderFactory(real sizeFactor, unsigned numFaces);
+	~HeuristicGridDensityProviderFactory ();
+
+	auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source, const real proscenium[4]);
+	auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform);
+	auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source);
+
+protected:
+	real sizeFactor;
+	unsigned numFaces;
+};
+
+#endif // HEURISTICGRIDDENSITYPROVIDERFACTORY_H
+
diff --git a/source/blender/freestyle/intern/view_map/OccluderSource.cpp b/source/blender/freestyle/intern/view_map/OccluderSource.cpp
new file mode 100644
index 00000000000..356e281be4b
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/OccluderSource.cpp
@@ -0,0 +1,132 @@
+//
+//  Filename         : OccluderSource.h
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2010-12-21
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 "OccluderSource.h"
+#include <algorithm>
+
+OccluderSource::OccluderSource (const GridHelpers::Transform& t, WingedEdge& we) : wingedEdge(we), valid(false), transform(t) {
+	begin();
+}
+
+OccluderSource::~OccluderSource() {
+}
+
+void OccluderSource::buildCachedPolygon() {
+	vector<Vec3r> vertices(GridHelpers::enumerateVertices((*currentFace)->getEdgeList()));
+	// This doesn't work, because our functor's polymorphism won't survive the copy:
+	// std::transform(vertices.begin(), vertices.end(), vertices.begin(), transform);
+	// so we have to do:
+	for ( vector<Vec3r>::iterator i = vertices.begin(); i != vertices.end(); ++i ) {
+		(*i) = transform(*i);
+	}
+	cachedPolygon = Polygon3r(vertices, transform((*currentFace)->GetNormal()));
+}
+
+void OccluderSource::begin() {
+ 	vector<WShape*>& wshapes = wingedEdge.getWShapes();
+	currentShape = wshapes.begin();
+	shapesEnd = wshapes.end();
+	valid = false;
+	if ( currentShape != shapesEnd ) {
+		vector<WFace*>& wFaces = (*currentShape)->GetFaceList();
+		currentFace = wFaces.begin();
+		facesEnd = wFaces.end();
+
+		if ( currentFace != facesEnd ) {
+			buildCachedPolygon();
+			valid = true;
+		}
+	}
+}
+
+bool OccluderSource::next() {
+	if ( valid ) {
+		++currentFace;
+		while ( currentFace == facesEnd ) {
+			++currentShape;
+			if ( currentShape == shapesEnd ) {
+				valid = false;
+				return false;
+			} else {
+				vector<WFace*>& wFaces = (*currentShape)->GetFaceList();
+				currentFace = wFaces.begin();
+				facesEnd = wFaces.end();
+			}
+		}
+		buildCachedPolygon();
+		return true;
+	}
+	return false;
+}
+
+bool OccluderSource::isValid() {
+	// Or:
+	// return currentShapes != shapesEnd && currentFace != facesEnd;
+	return valid;
+}
+
+WFace* OccluderSource::getWFace() {
+	return valid ? *currentFace : NULL;
+}
+
+Polygon3r OccluderSource::getCameraSpacePolygon() {
+	return Polygon3r(GridHelpers::enumerateVertices((*currentFace)->getEdgeList()), (*currentFace)->GetNormal());
+}
+
+Polygon3r& OccluderSource::getGridSpacePolygon() {
+	return cachedPolygon;
+}
+
+void OccluderSource::getOccluderProscenium(real proscenium[4]) {
+	begin();
+	const Vec3r& initialPoint = cachedPolygon.getVertices()[0];
+	proscenium[0] = proscenium[1] = initialPoint[0];
+	proscenium[2] = proscenium[3] = initialPoint[1];
+	while ( isValid() ) {
+		GridHelpers::expandProscenium (proscenium, cachedPolygon);
+		next();
+	}
+	cout << "Proscenium: (" << proscenium[0] << ", " << proscenium[1] << ", " << proscenium[2] << ", " << proscenium[3] << ")" << endl;
+}
+
+real OccluderSource::averageOccluderArea() {
+	real area = 0.0;
+	unsigned numFaces = 0;
+	for ( begin(); isValid(); next() ) {
+		Vec3r min, max;
+		cachedPolygon.getBBox(min, max);
+		area += (max[0] - min[0]) * (max[1] - min[1]);
+		++numFaces;
+	}
+	area /= numFaces;
+	return area;
+}
+
+
diff --git a/source/blender/freestyle/intern/view_map/OccluderSource.h b/source/blender/freestyle/intern/view_map/OccluderSource.h
new file mode 100644
index 00000000000..d623cc63fed
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/OccluderSource.h
@@ -0,0 +1,70 @@
+//
+//  Filename         : OccluderSource.h
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2010-12-21
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 OCCLUDERSOURCE_H
+#define OCCLUDERSOURCE_H
+
+#include "../winged_edge/WEdge.h"
+#include "../geometry/GridHelpers.h"
+
+class OccluderSource {
+	// Disallow copying and assignment
+	OccluderSource (const OccluderSource& other);
+	OccluderSource& operator= (const OccluderSource& other);
+
+public:
+	OccluderSource (const GridHelpers::Transform& transform, WingedEdge& we);
+	virtual ~OccluderSource();
+
+	void begin();
+	virtual bool next();
+	bool isValid();
+
+	WFace* getWFace();
+	Polygon3r getCameraSpacePolygon();
+	Polygon3r& getGridSpacePolygon();
+
+	virtual void getOccluderProscenium(real proscenium[4]);
+	virtual real averageOccluderArea();
+
+protected:
+	WingedEdge& wingedEdge;
+	vector<WShape*>::const_iterator currentShape, shapesEnd;
+	vector<WFace*>::const_iterator currentFace, facesEnd;
+
+	bool valid;
+
+	Polygon3r cachedPolygon;
+	const GridHelpers::Transform& transform;
+
+	void buildCachedPolygon();
+};
+
+#endif // OCCLUDERSOURCE_H
diff --git a/source/blender/freestyle/intern/view_map/Pow23GridDensityProvider.cpp b/source/blender/freestyle/intern/view_map/Pow23GridDensityProvider.cpp
new file mode 100644
index 00000000000..2e506da9ee9
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/Pow23GridDensityProvider.cpp
@@ -0,0 +1,109 @@
+//
+//  Filename         : Pow23GridDensityProvider.cpp
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2011-2-8
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 "Pow23GridDensityProvider.h"
+
+Pow23GridDensityProvider::Pow23GridDensityProvider(OccluderSource& source, const real proscenium[4], unsigned numFaces) 
+	: GridDensityProvider(source), numFaces(numFaces)
+{
+	initialize (proscenium);
+}
+
+Pow23GridDensityProvider::Pow23GridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform, unsigned numFaces) 
+	: GridDensityProvider(source), numFaces(numFaces)
+{
+	real proscenium[4];
+	calculateQuickProscenium(transform, bbox, proscenium);
+	
+	initialize (proscenium);
+}
+
+Pow23GridDensityProvider::Pow23GridDensityProvider(OccluderSource& source, unsigned numFaces) 
+	: GridDensityProvider(source), numFaces(numFaces) 
+{
+	real proscenium[4];
+	calculateOptimalProscenium(source, proscenium);
+
+	initialize (proscenium);
+}
+
+Pow23GridDensityProvider::~Pow23GridDensityProvider () {}
+
+void Pow23GridDensityProvider::initialize (const real proscenium[4]) 
+{
+	float prosceniumWidth = (proscenium[1] - proscenium[0]);
+	float prosceniumHeight = (proscenium[3] - proscenium[2]);
+	real cellArea = prosceniumWidth * prosceniumHeight / pow(numFaces, 2.0f / 3.0f);
+	cout << prosceniumWidth << " x " << prosceniumHeight << " grid with cells of area " << cellArea << "." << endl;
+
+	_cellSize = sqrt(cellArea);
+	// Now we know how many cells make each side of our grid
+	_cellsX = ceil(prosceniumWidth / _cellSize);
+	_cellsY = ceil(prosceniumHeight / _cellSize);
+	cout << _cellsX << "x" << _cellsY << " cells of size " << _cellSize << " square." << endl;
+
+	// Make sure the grid exceeds the proscenium by a small amount
+	float safetyZone = 0.1;
+	if ( _cellsX * _cellSize < prosceniumWidth * (1.0 + safetyZone) ) {
+		_cellsX = prosceniumWidth * (1.0 + safetyZone) / _cellSize;
+	}
+	if ( _cellsY * _cellSize < prosceniumHeight * (1.0 + safetyZone) ) {
+		_cellsY = prosceniumHeight * (1.0 + safetyZone) / _cellSize;
+	}
+	cout << _cellsX << "x" << _cellsY << " cells of size " << _cellSize << " square." << endl;
+
+	// Find grid origin
+	_cellOrigin[0] = ((proscenium[0] + proscenium[1]) / 2.0) - (_cellsX / 2.0) * _cellSize;
+	_cellOrigin[1] = ((proscenium[2] + proscenium[3]) / 2.0) - (_cellsY / 2.0) * _cellSize;
+}
+
+Pow23GridDensityProviderFactory::Pow23GridDensityProviderFactory(unsigned numFaces)
+	: numFaces(numFaces)
+{
+}
+
+Pow23GridDensityProviderFactory::~Pow23GridDensityProviderFactory () {}
+
+auto_ptr<GridDensityProvider> Pow23GridDensityProviderFactory::newGridDensityProvider(OccluderSource& source, const real proscenium[4]) 
+{
+	return auto_ptr<GridDensityProvider>(new Pow23GridDensityProvider(source, proscenium, numFaces));
+}
+
+auto_ptr<GridDensityProvider> Pow23GridDensityProviderFactory::newGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform) 
+{
+	return auto_ptr<GridDensityProvider>(new Pow23GridDensityProvider(source, bbox, transform, numFaces));
+}
+
+auto_ptr<GridDensityProvider> Pow23GridDensityProviderFactory::newGridDensityProvider(OccluderSource& source) 
+{
+	return auto_ptr<GridDensityProvider>(new Pow23GridDensityProvider(source, numFaces));
+}
+
+
diff --git a/source/blender/freestyle/intern/view_map/Pow23GridDensityProvider.h b/source/blender/freestyle/intern/view_map/Pow23GridDensityProvider.h
new file mode 100644
index 00000000000..c398a1643a4
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/Pow23GridDensityProvider.h
@@ -0,0 +1,67 @@
+//
+//  Filename         : Pow23GridDensityProvider.h
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2011-2-8
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 POW23GRIDDENSITYPROVIDER_H
+#define POW23GRIDDENSITYPROVIDER_H
+
+#include "GridDensityProvider.h"
+
+class Pow23GridDensityProvider : public GridDensityProvider {
+	// Disallow copying and assignment
+	Pow23GridDensityProvider (const Pow23GridDensityProvider& other);
+	Pow23GridDensityProvider& operator= (const Pow23GridDensityProvider& other);
+
+public:
+	Pow23GridDensityProvider(OccluderSource& source, const real proscenium[4], unsigned numFaces);
+	Pow23GridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform, unsigned numFaces);
+	Pow23GridDensityProvider(OccluderSource& source, unsigned numFaces);
+	virtual ~Pow23GridDensityProvider ();
+
+protected:
+	unsigned numFaces;
+
+private:
+	void initialize (const real proscenium[4]); 
+};
+
+class Pow23GridDensityProviderFactory : public GridDensityProviderFactory {
+public:
+	Pow23GridDensityProviderFactory(unsigned numFaces);
+	~Pow23GridDensityProviderFactory ();
+
+	auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source, const real proscenium[4]);
+	auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source, const BBox<Vec3r>& bbox, const GridHelpers::Transform& transform);
+	auto_ptr<GridDensityProvider> newGridDensityProvider(OccluderSource& source);
+protected:
+	unsigned numFaces;
+};
+
+#endif // POW23GRIDDENSITYPROVIDER_H
+
diff --git a/source/blender/freestyle/intern/view_map/Silhouette.h b/source/blender/freestyle/intern/view_map/Silhouette.h
index 86a4efcdc48..da337c1d0b3 100755
--- a/source/blender/freestyle/intern/view_map/Silhouette.h
+++ b/source/blender/freestyle/intern/view_map/Silhouette.h
@@ -404,6 +404,8 @@ protected:
 
   bool _isSmooth;
 
+  bool _isInImage;
+
 public:
   /*! A field that can be used by the user to store any data.
    *  This field must be reseted afterwards using ResetUserData().
@@ -421,6 +423,7 @@ public:
     //_hasVisibilityPoint=false;
     _occludeeEmpty = true;
     _isSmooth = false;
+	_isInImage = true;
   }
   /*! Builds an FEdge going from vA to vB. */  
   inline FEdge(SVertex *vA, SVertex *vB) {
@@ -434,6 +437,7 @@ public:
     //_hasVisibilityPoint=false;
     _occludeeEmpty = true;
     _isSmooth = false;
+	_isInImage = true;
   }
   /*! Copy constructor */
   inline FEdge(FEdge& iBrother)
@@ -451,6 +455,7 @@ public:
     _aFace = iBrother._aFace;
     _occludeeEmpty = iBrother._occludeeEmpty;
     _isSmooth = iBrother._isSmooth;
+	_isInImage = iBrother._isInImage;
     iBrother.userdata = this;
     userdata = 0;
   }
@@ -498,6 +503,7 @@ public:
   inline bool getOccludeeEmpty() { return _occludeeEmpty; }
   /*! Returns true if this FEdge is a smooth FEdge. */
   inline bool isSmooth() const {return _isSmooth;}
+  inline bool isInImage () const { return _isInImage; }
 
   /* modifiers */
   /*! Sets the first SVertex. */
@@ -525,6 +531,7 @@ public:
    *  true for Smooth, false for Sharp.
    */
   inline void setSmooth(bool iFlag) {_isSmooth = iFlag;}
+  inline void setIsInImage (bool iFlag) { _isInImage = iFlag; }
   
   /* checks whether two FEdge have a common vertex.
    *  Returns a pointer on the common vertex if it exists, 
diff --git a/source/blender/freestyle/intern/view_map/SphericalGrid.cpp b/source/blender/freestyle/intern/view_map/SphericalGrid.cpp
new file mode 100644
index 00000000000..c1647050d0d
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/SphericalGrid.cpp
@@ -0,0 +1,221 @@
+//
+//  Filename         : SphericalGrid.h
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2010-12-19
+//
+///////////////////////////////////////////////////////////////////////////////
+
+//
+//  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 "SphericalGrid.h"
+
+#include <stdexcept>
+#include <algorithm>
+
+using namespace std;
+
+// Helper Classes
+
+// OccluderData
+///////////////
+
+// Cell
+/////////
+
+SphericalGrid::Cell::Cell () {}
+
+SphericalGrid::Cell::~Cell () {}
+
+void SphericalGrid::Cell::setDimensions(real x, real y, real sizeX, real sizeY) {
+	const real epsilon = 1.0e-06;
+	boundary[0] = x - epsilon;
+	boundary[1] = x + sizeX + epsilon;
+	boundary[2] = y - epsilon;
+	boundary[3] = y + sizeY + epsilon;
+}
+
+bool SphericalGrid::Cell::compareOccludersByShallowestPoint (const SphericalGrid::OccluderData* a, const SphericalGrid::OccluderData* b) {
+	return a->shallowest < b->shallowest;
+}
+
+void SphericalGrid::Cell::indexPolygons() {
+	// Sort occluders by their shallowest points.
+	sort(faces.begin(), faces.end(), compareOccludersByShallowestPoint);
+}
+
+// Iterator
+//////////////////
+
+SphericalGrid::Iterator::Iterator (SphericalGrid& grid, Vec3r& center, real epsilon) 
+	: _target(SphericalGrid::Transform::sphericalProjection(center)), 
+	_foundOccludee(false)
+{
+	// Find target cell
+	_cell = grid.findCell(_target);
+	#if sphericalgridlogging == 1
+		cout << "Searching for occluders of edge centered at " << _target << " in cell [" 
+			<< _cell->boundary[0] << ", " << _cell->boundary[1] << ", " << _cell->boundary[2] 
+			<< ", " << _cell->boundary[3] << "] (" << _cell->faces.size() << " occluders)" << endl;
+	#endif
+
+	// Set iterator
+	_current = _cell->faces.begin();
+}
+
+SphericalGrid::Iterator::~Iterator () {}
+
+// SphericalGrid
+/////////////////
+
+SphericalGrid::SphericalGrid(OccluderSource& source, GridDensityProvider& density, ViewMap *viewMap, Vec3r& viewpoint, bool enableQI)
+	: _viewpoint(viewpoint),
+	_enableQI(enableQI)
+{
+	cout << "Generate Cell structure" << endl;
+	// Generate Cell structure
+	assignCells(source, density, viewMap);
+	cout << "Distribute occluders" << endl;
+	// Fill Cells
+	distributePolygons(source);
+	cout << "Reorganize cells" << endl;
+	// Reorganize Cells
+	reorganizeCells();
+	cout << "Ready to use SphericalGrid" << endl;
+}
+
+SphericalGrid::~SphericalGrid () {
+}
+
+void SphericalGrid::assignCells (OccluderSource& source, GridDensityProvider& density, ViewMap *viewMap) {
+	_cellSize = density.cellSize();
+	_cellsX = density.cellsX();
+	_cellsY = density.cellsY();
+	_cellOrigin[0] = density.cellOrigin(0);
+	_cellOrigin[1] = density.cellOrigin(1);
+
+	// Now allocate the cell table and fill it with default (empty) cells
+	_cells.resize(_cellsX * _cellsY);
+	for ( cellContainer::iterator i = _cells.begin(), end = _cells.end(); i != end; ++i ) {
+		(*i) = NULL;
+	}
+
+	// Identify cells that will be used, and set the dimensions for each
+	ViewMap::fedges_container& fedges = viewMap->FEdges();
+	for (ViewMap::fedges_container::iterator f = fedges.begin(), fend = fedges.end(); f != fend; ++f ) {
+		if ( (*f)->isInImage() ) {
+			Vec3r point = SphericalGrid::Transform::sphericalProjection((*f)->center3d());
+			unsigned i, j;
+			getCellCoordinates(point, i, j);
+			if ( _cells[i * _cellsY + j] == NULL ) {
+				// This is an uninitialized cell
+
+				real x, y, width, height;
+
+				x = _cellOrigin[0] + _cellSize * i;
+				width = _cellSize;
+
+				y = _cellOrigin[1] + _cellSize * j;
+				height = _cellSize;
+
+				// Initialize cell
+				Cell* b = _cells[i * _cellsY + j] = new Cell();
+				b->setDimensions(x, y, width, height);
+			}
+		}
+	}
+}
+
+void SphericalGrid::distributePolygons (OccluderSource& source) {
+	unsigned long nFaces = 0;
+	unsigned long nKeptFaces = 0;
+
+	for ( source.begin(); source.isValid(); source.next() ) {
+		OccluderData* occluder = NULL;
+
+		try {
+			if ( insertOccluder(source, occluder) ) {
+				_faces.push_back(occluder);
+				++nKeptFaces;
+			}
+		} catch (...) {
+			// If an exception was thrown, _faces.push_back() cannot have succeeded.
+			// occluder is not owned by anyone, and must be deleted.
+			// If the exception was thrown before or during new OccluderData(), then
+			// occluder is NULL, and this delete is harmless.
+			delete occluder;
+			throw;
+		}
+		++nFaces;
+	}
+	cout << "Distributed " << nFaces << " occluders.  Retained " << nKeptFaces << "." << endl;
+}
+
+void SphericalGrid::reorganizeCells () {
+	// Sort the occluders by shallowest point
+	for ( vector<Cell*>::iterator i = _cells.begin(), end = _cells.end(); i != end; ++i ) {
+		if ( *i != NULL ) {
+			(*i)->indexPolygons();
+		}
+	}
+}
+
+void SphericalGrid::getCellCoordinates(const Vec3r& point, unsigned& x, unsigned& y) {
+	x = min(_cellsX - 1, (unsigned) floor (max((double) 0.0f, point[0] - _cellOrigin[0]) / _cellSize));
+	y = min(_cellsY - 1, (unsigned) floor (max((double) 0.0f, point[1] - _cellOrigin[1]) / _cellSize));
+}
+
+SphericalGrid::Cell* SphericalGrid::findCell(const Vec3r& point) {
+	unsigned x, y;
+	getCellCoordinates(point, x, y);
+	return _cells[x * _cellsY + y];
+}
+
+bool SphericalGrid::orthographicProjection () const {
+	return false;
+}
+
+const Vec3r& SphericalGrid::viewpoint() const {
+	return _viewpoint;
+}
+
+bool SphericalGrid::enableQI() const {
+	return _enableQI;
+}
+
+SphericalGrid::Transform::Transform () : GridHelpers::Transform() {
+}
+
+Vec3r SphericalGrid::Transform::operator() (const Vec3r& point) const {
+	return sphericalProjection(point);
+}
+
+Vec3r SphericalGrid::Transform::sphericalProjection(const Vec3r& M) {
+	Vec3r newPoint;
+
+	newPoint[0] = ::atan(M[0] / M[2]);
+	newPoint[1] = ::atan(M[1] / M[2]);
+	newPoint[2] = ::sqrt(M[0] * M[0] + M[1] * M[1] + M[2] * M[2]);
+
+	return newPoint;
+}
+
diff --git a/source/blender/freestyle/intern/view_map/SphericalGrid.h b/source/blender/freestyle/intern/view_map/SphericalGrid.h
new file mode 100644
index 00000000000..55cdee19a80
--- /dev/null
+++ b/source/blender/freestyle/intern/view_map/SphericalGrid.h
@@ -0,0 +1,388 @@
+//
+//  Filename         : SphericalGrid.h
+//  Author(s)        : Alexander Beels
+//  Purpose          : Class to define a cell grid surrounding
+//                     the projected image of a scene
+//  Date of creation : 2010-12-19
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+//
+//  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 SPHERICALGRID_H
+#define SPHERICALGRID_H
+
+#define sphericalgridlogging 0
+
+// I would like to avoid using deque because including ViewMap.h and <deque> or <vector>
+// separately results in redefinitions of identifiers.  ViewMap.h already includes <vector> 
+// so it should be a safe fall-back.
+//#include <vector>
+//#include <deque>
+#include "ViewMap.h"
+#include "../winged_edge/WEdge.h"
+#include "../geometry/Polygon.h"
+#include "../system/PointerSequence.h"
+#include "../geometry/BBox.h"
+#include "../geometry/GridHelpers.h"
+#include "OccluderSource.h"
+#include "GridDensityProvider.h"
+
+class SphericalGrid
+{
+public:
+	// Helper classes
+	struct OccluderData {
+		explicit OccluderData (OccluderSource& source, Polygon3r& p);
+		Polygon3r poly;
+		Polygon3r cameraSpacePolygon;
+		real shallowest, deepest;
+		// N.B. We could, of course, store face in poly's userdata
+		// member, like the old ViewMapBuilder code does.  However,
+		// code comments make it clear that userdata is deprecated,
+		// so we avoid the temptation to save 4 or 8 bytes.
+		WFace* face;
+	};
+
+private:
+	struct Cell {
+		// Can't store Cell in a vector without copy and assign
+		//Cell(const Cell& other);
+		//Cell& operator= (const Cell& other);
+
+		explicit Cell ();
+		~Cell ();
+
+		static bool compareOccludersByShallowestPoint (const OccluderData* a, const OccluderData* b);
+
+		void setDimensions(real x, real y, real sizeX, real sizeY);
+		void checkAndInsert(OccluderSource& source, Polygon3r& poly, OccluderData*& occluder);
+		void indexPolygons();
+
+		real boundary[4];
+		//deque<OccluderData*> faces;
+		vector<OccluderData*> faces;
+	};
+
+public:
+	/*****
+
+	Iterator needs to allow the user to avoid full 3D comparison in
+	two cases:
+
+	(1) Where (*current)->deepest < target[2], where the occluder is 
+	unambiguously in front of the target point.
+
+	(2) Where (*current)->shallowest > target[2], where the occluder
+	is unambiguously in back of the target point.
+
+	In addition, when used by OptimizedFindOccludee, Iterator should
+	stop iterating as soon as it has an occludee candidate and 
+	(*current)->shallowest > candidate[2], because at that point forward
+	no new occluder could possibly be a better occludee.
+
+	*****/
+
+	class Iterator {
+	public:
+		// epsilon is not used in this class, but other grids with the same interface may need an epsilon
+		explicit Iterator (SphericalGrid& grid, Vec3r& center, real epsilon=1e-06);
+		~Iterator ();
+		void initBeforeTarget ();
+		void initAfterTarget ();
+		void nextOccluder ();
+		void nextOccludee ();
+		bool validBeforeTarget();
+		bool validAfterTarget();
+		WFace* getWFace() const;
+		Polygon3r* getCameraSpacePolygon();
+		void reportDepth(Vec3r origin, Vec3r u, real t);
+	private:
+		bool testOccluder(bool wantOccludee);
+		void markCurrentOccludeeCandidate(real depth);
+
+		Cell* _cell;
+		Vec3r _target;
+		bool _foundOccludee;
+		real _occludeeDepth;
+		//deque<OccluderData*>::iterator _current, _occludeeCandidate;
+		vector<OccluderData*>::iterator _current, _occludeeCandidate;
+	};
+
+	class Transform : public GridHelpers::Transform {
+	public:
+		explicit Transform ();
+		explicit Transform (Transform& other);
+		Vec3r operator() (const Vec3r& point) const;
+		static Vec3r sphericalProjection(const Vec3r& M);
+	};
+
+private:
+	// Prevent implicit copies and assignments.
+	SphericalGrid(const SphericalGrid& other);
+	SphericalGrid& operator= (const SphericalGrid& other);
+public:
+	explicit SphericalGrid (OccluderSource& source, GridDensityProvider& density, ViewMap *viewMap, Vec3r& viewpoint, bool enableQI);
+	virtual ~SphericalGrid();
+
+	// Generate Cell structure
+	void assignCells(OccluderSource& source, GridDensityProvider& density, ViewMap *viewMap);
+	// Fill Cells
+	void distributePolygons(OccluderSource& source);
+	// Insert one polygon into each matching cell, 
+	// return true if any cell consumes the polygon
+	bool insertOccluder(OccluderSource& source, OccluderData*& occluder);
+	// Sort occluders in each cell
+	void reorganizeCells();
+
+	Cell* findCell(const Vec3r& point);
+
+	// Accessors:
+	bool orthographicProjection() const;
+	const Vec3r& viewpoint() const;
+	bool enableQI() const; 
+
+private:
+	void getCellCoordinates(const Vec3r& point, unsigned& x, unsigned& y);
+
+	typedef PointerSequence<vector<Cell*>, Cell*> cellContainer;
+	//typedef PointerSequence<deque<OccluderData*>, OccluderData*> occluderContainer;
+	typedef PointerSequence<vector<OccluderData*>, OccluderData*> occluderContainer;
+	unsigned _cellsX, _cellsY;
+	float _cellSize;
+	float _cellOrigin[2];
+	cellContainer _cells;
+	occluderContainer _faces;
+	Vec3r _viewpoint;
+	bool _enableQI;
+};
+
+inline void SphericalGrid::Iterator::initBeforeTarget () {
+	_current = _cell->faces.begin();
+	while ( _current != _cell->faces.end() && ! testOccluder(false) ) {
+		++_current;
+	}
+}
+
+inline void SphericalGrid::Iterator::initAfterTarget () {
+	if ( _foundOccludee ) {
+		#if sphericalgridlogging == 1
+			std::cout << "\tStarting occludee search from occludeeCandidate at depth " << _occludeeDepth << std::endl;
+		#endif
+		_current = _occludeeCandidate;
+		return;
+	}
+
+	#if sphericalgridlogging == 1
+		std::cout << "\tStarting occludee search from current position" << std::endl;
+	#endif
+
+	while ( _current != _cell->faces.end() && ! testOccluder(true) ) {
+		++_current;
+	}
+}
+
+inline bool SphericalGrid::Iterator::testOccluder (bool wantOccludee) {
+	// End-of-list is not even a valid iterator position
+	if ( _current == _cell->faces.end() ) {
+		// Returning true seems strange, but it will break us out of whatever loop
+		// is calling testOccluder, and _current=_cell->face.end() will make
+		// the calling routine give up.
+		return true;
+	}
+	#if sphericalgridlogging == 1
+		std::cout << "\tTesting occluder " << (*_current)->poly.getVertices()[0];
+		for ( unsigned i = 1; i < (*_current)->poly.getVertices().size(); ++i ) {
+			std::cout << ", " << (*_current)->poly.getVertices()[i];
+		}
+		std::cout << " from shape " << (*_current)->face->GetVertex(0)->shape()->GetId() << std::endl;
+	#endif
+
+
+	// If we have an occluder candidate and we are unambiguously after it, abort
+	if ( _foundOccludee && (*_current)->shallowest > _occludeeDepth ) {
+		#if sphericalgridlogging == 1
+			std::cout << "\t\tAborting: shallowest > occludeeCandidate->deepest" << std::endl;
+		#endif
+		_current = _cell->faces.end();
+
+		// See note above
+		return true;
+	}
+
+	// Specific continue or stop conditions when searching for each type
+	if ( wantOccludee ) {
+		if ( (*_current)->deepest < _target[2] ) {
+			#if sphericalgridlogging == 1
+				std::cout << "\t\tSkipping: shallower than target while looking for occludee" << std::endl;
+			#endif
+			return false;
+		}
+	} else {
+		if ( (*_current)->shallowest > _target[2] ) {
+			#if sphericalgridlogging == 1
+				std::cout << "\t\tStopping: deeper than target while looking for occluder" << std::endl;
+			#endif
+			return true;
+		}
+	}
+
+	// Depthwise, this is a valid occluder.
+
+	// Check to see if target is in the 2D bounding box
+	Vec3r bbMin, bbMax;
+	(*_current)->poly.getBBox(bbMin, bbMax);
+	if ( _target[0] < bbMin[0] || _target[0] > bbMax[0] || _target[1] < bbMin[1] || _target[1] > bbMax[1] ) {
+		#if sphericalgridlogging == 1
+			std::cout << "\t\tSkipping: bounding box violation" << std::endl;
+		#endif
+		return false;
+	}
+
+	// We've done all the corner cutting we can.
+	// Let the caller work out whether or not
+	// the geometry is correct.
+	return true;
+}
+
+inline void SphericalGrid::Iterator::reportDepth (Vec3r origin, Vec3r u, real t) {
+	// The reported depth is the length of a ray in camera space
+	// We need to convert it into the distance from viewpoint
+	// If origin is the viewpoint, depth == t
+	// A future optimization could allow the caller to tell us if origin is viewponit or target,
+	// at the cost of changing the OptimizedGrid API.
+	real depth = (origin + u * t).norm();
+	#if sphericalgridlogging == 1
+		std::cout << "\t\tReporting depth of occluder/ee: " << depth;
+	#endif
+	if ( depth > _target[2] ) {
+		#if sphericalgridlogging == 1
+			std::cout << " is deeper than target" << std::endl;
+		#endif
+		// If the current occluder is the best occludee so far, save it.
+		if ( ! _foundOccludee || _occludeeDepth > depth ) {
+			markCurrentOccludeeCandidate(depth);
+		} 
+	} else {
+		#if sphericalgridlogging == 1
+			std::cout << std::endl;
+		#endif
+	}
+}
+
+inline void SphericalGrid::Iterator::nextOccluder () {
+	if ( _current != _cell->faces.end() ) {
+		do {
+			++_current;
+		} while ( _current != _cell->faces.end() && ! testOccluder(false) );
+	}
+}
+
+inline void SphericalGrid::Iterator::nextOccludee () {
+	if ( _current != _cell->faces.end() ) {
+		do {
+			++_current;
+		} while ( _current != _cell->faces.end() && ! testOccluder(true) );
+	}
+}
+
+inline bool SphericalGrid::Iterator::validBeforeTarget () {
+	return _current != _cell->faces.end() && (*_current)->shallowest <= _target[2];
+}
+
+inline bool SphericalGrid::Iterator::validAfterTarget () {
+	return _current != _cell->faces.end();
+}
+
+inline void SphericalGrid::Iterator::markCurrentOccludeeCandidate(real depth) {
+	#if sphericalgridlogging == 1
+		std::cout << "\t\tFound occludeeCandidate at depth " << depth << std::endl;
+	#endif
+	_occludeeCandidate = _current;
+	_occludeeDepth = depth;
+	_foundOccludee = true;
+}
+
+inline WFace* SphericalGrid::Iterator::getWFace() const {
+	return (*_current)->face;
+}
+
+inline Polygon3r* SphericalGrid::Iterator::getCameraSpacePolygon() {
+	return &((*_current)->cameraSpacePolygon);
+}
+
+inline SphericalGrid::OccluderData::OccluderData (OccluderSource& source, Polygon3r& p)
+	: poly(p), 
+	cameraSpacePolygon(source.getCameraSpacePolygon()),
+	face(source.getWFace())
+{
+	const Vec3r viewpoint(0, 0, 0);
+	// Get the point on the camera-space polygon that is closest to the viewpoint
+	// shallowest is the distance from the viewpoint to that point
+	shallowest = GridHelpers::distancePointToPolygon(viewpoint, cameraSpacePolygon);
+
+	// Get the point on the camera-space polygon that is furthest from the viewpoint
+	// deepest is the distance from the viewpoint to that point
+	deepest = cameraSpacePolygon.getVertices()[2].norm();
+	for ( unsigned i = 0; i < 2; ++i ) {
+		real t = cameraSpacePolygon.getVertices()[i].norm();
+		if ( t > deepest ) {
+			deepest = t;
+		}
+	}
+}
+
+inline void SphericalGrid::Cell::checkAndInsert(OccluderSource& source, Polygon3r& poly, OccluderData*& occluder) {
+	if ( GridHelpers::insideProscenium (boundary, poly) ) {
+		if ( occluder == NULL) {
+			// Disposal of occluder will be handled in SphericalGrid::distributePolygons(),
+			// or automatically by SphericalGrid::_faces;
+			occluder = new OccluderData(source, poly);
+		}
+		faces.push_back(occluder);
+	}
+}
+
+inline bool SphericalGrid::insertOccluder(OccluderSource& source, OccluderData*& occluder) {
+	Polygon3r& poly(source.getGridSpacePolygon());
+	occluder = NULL;
+
+	Vec3r bbMin, bbMax;
+	poly.getBBox(bbMin, bbMax);
+	// Check overlapping cells
+	unsigned startX, startY, endX, endY;
+	getCellCoordinates(bbMin, startX, startY);
+	getCellCoordinates(bbMax, endX, endY);
+
+	for ( unsigned i = startX; i <= endX; ++i ) {
+		for ( unsigned j = startY; j <= endY; ++j ) {
+			if ( _cells[i * _cellsY + j] != NULL ) {
+				_cells[i * _cellsY + j]->checkAndInsert(source, poly, occluder);
+			}
+		}
+	}
+
+	return occluder != NULL;
+}
+
+#endif // SPHERICALGRID_H
+
diff --git a/source/blender/freestyle/intern/view_map/ViewMap.h b/source/blender/freestyle/intern/view_map/ViewMap.h
index a0e60ff7bd7..d2d7a9a5341 100755
--- a/source/blender/freestyle/intern/view_map/ViewMap.h
+++ b/source/blender/freestyle/intern/view_map/ViewMap.h
@@ -801,6 +801,7 @@ private:
   // and _aShape is the one on its right // NON GERE PAR LE COPY CONSTRUCTEUR
   int _qi;
   vector<ViewShape*> _Occluders;
+  bool _isInImage;
 
   // tmp
   Id * _splittingId;
@@ -821,6 +822,7 @@ public:
     _aShape=0;
     userdata = 0;
     _splittingId = 0;
+    _isInImage = true;
   }
   inline ViewEdge(ViewVertex* iA, ViewVertex *iB)
   {
@@ -834,6 +836,7 @@ public:
     _qi = 0;
     userdata = 0;
     _splittingId = 0;
+    _isInImage = true;
   }
   inline ViewEdge(ViewVertex* iA, ViewVertex *iB, FEdge *iFEdgeA)
   {
@@ -847,6 +850,7 @@ public:
     _qi = 0;
     userdata = 0;
     _splittingId = 0;
+    _isInImage = true;
   }
   inline ViewEdge(ViewVertex* iA, ViewVertex *iB, FEdge *iFEdgeA, FEdge *iFEdgeB, ViewShape *iShape)
   {
@@ -860,6 +864,7 @@ public:
     _qi = 0;
     userdata = 0;
     _splittingId = 0;
+    _isInImage = true;
     UpdateFEdges(); // tells every FEdge between iFEdgeA and iFEdgeB that this is theit ViewEdge
   }
 
@@ -879,6 +884,7 @@ public:
     _aShape = iBrother._aShape;
     _qi = iBrother._qi;
     _splittingId = 0;
+    _isInImage = iBrother._isInImage;
     iBrother.userdata = this;
     userdata = 0;
   }
@@ -937,6 +943,7 @@ public:
   inline const ViewShape * bShape() const {return _Shape;}
   inline vector<ViewShape*>& occluders() {return _Occluders;}
   inline Id * splittingId() {return _splittingId;}
+  inline bool isInImage() const { return _isInImage; }
   
   /* modifiers */
   /*! Sets the first ViewVertex of the ViewEdge. */
@@ -966,6 +973,7 @@ public:
   inline void setChainingTimeStamp(unsigned ts) {_ChainingTimeStamp = ts;}
   inline void AddOccluder(ViewShape *iShape) {_Occluders.push_back(iShape);}
   inline void setSplittingId(Id * id) {_splittingId = id;}
+  inline void setIsInImage(bool iFlag) { _isInImage = iFlag; }
 
   /* stroke interface definition */
   inline bool intersect_2d_area(const Vec2r& iMin, const Vec2r& iMax) const
diff --git a/source/blender/freestyle/intern/view_map/ViewMapBuilder.cpp b/source/blender/freestyle/intern/view_map/ViewMapBuilder.cpp
index dbd046f050a..242b769e8b5 100755
--- a/source/blender/freestyle/intern/view_map/ViewMapBuilder.cpp
+++ b/source/blender/freestyle/intern/view_map/ViewMapBuilder.cpp
@@ -21,15 +21,897 @@
 
 #include "ViewMapBuilder.h"
 #include <algorithm>
+#include <stdexcept>
+#include <memory>
+#include "../winged_edge/WFillGrid.h"
+#include "../../FRS_freestyle.h"
+#include "../geometry/GeomUtils.h"
+#include "../geometry/GridHelpers.h"
+#include "BoxGrid.h"
+#include "SphericalGrid.h"
+#include "OccluderSource.h"
+#include "CulledOccluderSource.h"
+#include "HeuristicGridDensityProviderFactory.h"
+
+#define logging 0
 
 using namespace std;
 
-ViewMap* ViewMapBuilder::BuildViewMap(WingedEdge& we, visibility_algo iAlgo, real epsilon) {
+template <typename G, typename I>
+static void findOccludee(FEdge *fe, G& grid, I& occluders, real epsilon, WFace** oaWFace, 
+    Vec3r& u, Vec3r& A, Vec3r& origin, Vec3r& edge, vector<WVertex*>& faceVertices)
+{
+  WFace *face = 0;
+  if(fe->isSmooth()){
+    FEdgeSmooth * fes = dynamic_cast<FEdgeSmooth*>(fe);
+    face = (WFace*)fes->face();
+  }
+  WFace * oface;
+  bool skipFace;
+  
+  WVertex::incoming_edge_iterator ie;
+  
+  *oaWFace = 0;
+  if(((fe)->getNature() & Nature::SILHOUETTE) || ((fe)->getNature() & Nature::BORDER))
+  {
+    // we cast a ray from A in the same direction but looking behind
+    Vec3r v(-u[0],-u[1],-u[2]);
+    bool noIntersection = true;
+    real mint=FLT_MAX;
+
+    for( occluders.initAfterTarget(); occluders.validAfterTarget(); occluders.nextOccludee() )
+    { 
+#if logging > 0
+	cout << "\t\tEvaluating intersection for occludee " << occluders.getWFace() << " and ray " << A << " * " << u << endl;
+#endif
+      oface = occluders.getWFace();
+      Polygon3r* p = occluders.getCameraSpacePolygon();
+      real d = -((p->getVertices())[0] * p->getNormal());
+      real t,t_u,t_v;
+      
+      if(0 != face)
+      { 
+        skipFace = false;
+        
+        if(face == oface)
+          continue;
+        
+        if(faceVertices.empty())
+          continue;
+        
+        for(vector<WVertex*>::iterator fv=faceVertices.begin(), fvend=faceVertices.end();
+        fv!=fvend;
+        ++fv)
+        { 
+          if((*fv)->isBoundary())
+            continue;
+          WVertex::incoming_edge_iterator iebegin=(*fv)->incoming_edges_begin();
+          WVertex::incoming_edge_iterator ieend=(*fv)->incoming_edges_end();
+          for(ie=iebegin;ie!=ieend; ++ie)
+          {  
+            if((*ie) == 0)
+              continue;
+            
+            WFace * sface = (*ie)->GetbFace();
+            if(sface == oface)
+            { 
+              skipFace = true;
+              break;
+            } 
+          }  
+          if(skipFace)
+            break;
+        } 
+        if(skipFace)
+          continue;
+      }  
+      else
+      {
+        // check whether the edge and the polygon plane are coincident:
+        //-------------------------------------------------------------
+        //first let us compute the plane equation.
+        if(GeomUtils::COINCIDENT == GeomUtils::intersectRayPlane(origin, edge, p->getNormal(), d, t, epsilon)) {
+#if logging > 0
+cout << "\t\tRejecting occluder for target coincidence." << endl;
+#endif
+          continue;
+        }
+      }
+
+      if(p->rayIntersect(A, v, t, t_u, t_v))
+      {
+#if logging > 0
+cout << "\t\tRay " << A << " * " << v << " intersects at time " << t << endl;
+#endif
+#if logging > 0
+cout << "\t\t(v * normal) == " << (v * p->getNormal()) << " for normal " << p->getNormal() << endl;
+#endif
+        if (fabs(v * p->getNormal()) > 0.0001)
+          if ((t>0.0)) // && (t<1.0))
+          {
+            if (t<mint)
+            {
+              *oaWFace = oface;
+              mint = t;
+              noIntersection = false;
+              fe->setOccludeeIntersection(Vec3r(A+t*v));
+#if logging > 0
+cout << "\t\tIs occludee" << endl;
+#endif
+            }
+          }
+
+        occluders.reportDepth(A, v, t);
+      }
+
+    }
+    
+    if(noIntersection)
+      *oaWFace = 0;
+  }
+}
+
+template <typename G, typename I>
+static void findOccludee(FEdge *fe, G& grid, real epsilon, ViewEdge* ve, WFace** oaFace)
+{
+  Vec3r A;
+  Vec3r edge;
+  Vec3r origin;
+  A = Vec3r(((fe)->vertexA()->point3D() + (fe)->vertexB()->point3D()) / 2.0);
+  edge = Vec3r((fe)->vertexB()->point3D()-(fe)->vertexA()->point3D());
+  origin = Vec3r((fe)->vertexA()->point3D());
+  Vec3r u;
+  if (grid.orthographicProjection()) {
+    u = Vec3r(0.0, 0.0, grid.viewpoint().z()-A.z());
+  } else {
+    u = Vec3r(grid.viewpoint()-A);
+  }
+  u.normalize();
+
+  vector<WVertex*> faceVertices;
+     
+  WFace *face = 0;
+  if(fe->isSmooth()) {
+    FEdgeSmooth * fes = dynamic_cast<FEdgeSmooth*>(fe);
+    face = (WFace*)fes->face();
+  }
+
+  if(0 != face) {
+    face->RetrieveVertexList(faceVertices);
+  }
+  
+  I occluders(grid, A, epsilon);
+  findOccludee<G, I>(fe, grid, occluders, epsilon, oaFace, u, A, origin, edge, faceVertices);
+}
+
+// computeVisibility takes a pointer to foundOccluders, instead of using a reference,
+// so that computeVeryFastVisibility can skip the AddOccluders step with minimal overhead.
+template <typename G, typename I>
+static int computeVisibility(ViewMap* viewMap, FEdge *fe, G& grid, real epsilon, ViewEdge* ve, WFace** oaWFace, set<ViewShape*>* foundOccluders)
+{
+  int qi = 0;
+
+  Vec3r center;
+  Vec3r edge;
+  Vec3r origin;
+
+  center = fe->center3d();
+  edge = Vec3r(fe->vertexB()->point3D() - fe->vertexA()->point3D());
+  origin = Vec3r(fe->vertexA()->point3D());
+
+  Vec3r vp;
+  if (grid.orthographicProjection()) {
+    vp = Vec3r(center.x(), center.y(), grid.viewpoint().z());
+  } else {
+    vp = Vec3r(grid.viewpoint());
+  }
+  Vec3r u(vp - center);
+  real raylength = u.norm();
+  u.normalize();
+  
+  WFace *face = 0;
+  if(fe->isSmooth()){
+    FEdgeSmooth * fes = dynamic_cast<FEdgeSmooth*>(fe);
+    face = (WFace*)fes->face();
+  }
+  vector<WVertex*> faceVertices;
+  WVertex::incoming_edge_iterator ie;
+            
+  WFace * oface;
+  bool skipFace;
+
+  if(face)
+    face->RetrieveVertexList(faceVertices);
+
+  I occluders(grid, center, epsilon);
+
+  for(occluders.initBeforeTarget(); occluders.validBeforeTarget(); occluders.nextOccluder())
+    { 
+      // If we're dealing with an exact silhouette, check whether 
+      // we must take care of this occluder of not.
+      // (Indeed, we don't consider the occluders that 
+      // share at least one vertex with the face containing 
+      // this edge).
+      //-----------
+      oface = occluders.getWFace();
+      Polygon3r* p = occluders.getCameraSpacePolygon();
+      real t, t_u, t_v;
+#if logging > 0
+	cout << "\t\tEvaluating intersection for occluder " << (p->getVertices())[0] << (p->getVertices())[1] << (p->getVertices())[2] << endl << "\t\t\tand ray " << vp << " * " << u << " (center " << center << ")" << endl;
+#endif
+
+#if logging > 0
+	Vec3r v(vp - center);
+	real rl = v.norm();
+	v.normalize();
+	vector<Vec3r> points;
+	// Iterate over vertices, storing projections in points
+	for(vector<WOEdge*>::const_iterator woe=oface->getEdgeList().begin(), woend=oface->getEdgeList().end(); woe!=woend; woe++) {
+		points.push_back(Vec3r((*woe)->GetaVertex()->GetVertex()));
+	}
+	Polygon3r p1(points, oface->GetNormal());
+	Vec3r v1((p1.getVertices())[0]);
+	real d = -(v1 * p->getNormal());
+	cout << "\t\tp:  " << (p->getVertices())[0] << (p->getVertices())[1] << (p->getVertices())[2] << ", norm: " << p->getNormal() << endl;
+	cout << "\t\tp1: " << (p1.getVertices())[0] << (p1.getVertices())[1] << (p1.getVertices())[2] << ", norm: " << p1.getNormal() << endl;
+#else
+      real d = -((p->getVertices())[0] * p->getNormal());
+#endif
+            
+      if(0 != face)
+	{
+#if logging > 0
+cout << "\t\tDetermining face adjacency...";
+#endif
+	  skipFace = false;
+              
+	  if(face == oface) {
+#if logging > 0
+cout << "  Rejecting occluder for face concurrency." << endl;
+#endif
+	    continue;
+          }
+              
+              
+	  for(vector<WVertex*>::iterator fv=faceVertices.begin(), fvend=faceVertices.end();
+              fv!=fvend;
+              ++fv)
+	    {
+                if((*fv)->isBoundary())
+                  continue;
+
+	      WVertex::incoming_edge_iterator iebegin=(*fv)->incoming_edges_begin();
+	      WVertex::incoming_edge_iterator ieend=(*fv)->incoming_edges_end();
+	      for(ie=iebegin;ie!=ieend; ++ie)
+                { 
+                  if((*ie) == 0)
+                    continue;
+                  
+                  WFace * sface = (*ie)->GetbFace();
+                  //WFace * sfacea = (*ie)->GetaFace();
+                  //if((sface == oface) || (sfacea == oface))
+                  if(sface == oface)
+		    {
+		      skipFace = true;
+		      break;
+		    }
+                }
+	      if(skipFace)
+		break;
+	    }
+	  if(skipFace) {
+#if logging > 0
+cout << "  Rejecting occluder for face adjacency." << endl;
+#endif
+	    continue;
+          }
+	}
+      else
+	{
+	  // check whether the edge and the polygon plane are coincident:
+	  //-------------------------------------------------------------
+	  //first let us compute the plane equation.
+           
+	  if(GeomUtils::COINCIDENT == GeomUtils::intersectRayPlane(origin, edge, p->getNormal(), d, t, epsilon)) {
+#if logging > 0
+cout << "\t\tRejecting occluder for target coincidence." << endl;
+#endif
+	    continue;
+          }
+	}
+
+#if logging > 0
+
+	real x;
+	if ( p1.rayIntersect(center, v, x, t_u, t_v) ) {
+		cout << "\t\tRay should intersect at time " << (rl - x) << ". Center: " << center << ", V: " << v <<  ", RL: " << rl << ", T:" << x << endl;
+	} else {
+		cout << "\t\tRay should not intersect. Center: " << center << ", V: " << v <<  ", RL: " << rl << endl;
+	}
+
+#endif
+
+      if(p->rayIntersect(center, u, t, t_u, t_v))
+	{
+#if logging > 0
+cout << "\t\tRay " << center << " * " << u << " intersects at time " << t << " (raylength is " << raylength << ")" << endl;
+#endif
+#if logging > 0
+cout << "\t\t(u * normal) == " << (u * p->getNormal()) << " for normal " << p->getNormal() << endl;
+#endif
+          if (fabs(u * p->getNormal()) > 0.0001)
+            if ((t>0.0) && (t<raylength))
+              {
+#if logging > 0
+cout << "\t\tIs occluder" << endl;
+#endif
+                if ( foundOccluders != NULL ) {
+                  ViewShape *vshape = viewMap->viewShape(oface->GetVertex(0)->shape()->GetId());
+                  foundOccluders->insert(vshape);
+                }
+
+                ++qi;
+
+                if(! grid.enableQI())
+                  break;
+              }
+
+          occluders.reportDepth(center, u, t);
+	}
+    }
+
+  // Find occludee
+  findOccludee<G, I>(fe, grid, occluders, epsilon, oaWFace, u, center, origin, edge, faceVertices);
+
+  return qi;
+}
+
+// computeCumulativeVisibility returns the lowest x such that the majority
+// of FEdges have QI <= x
+//
+// This was probably the original intention of the "normal" algorithm
+// on which computeDetailedVisibility is based.  But because the "normal"
+// algorithm chooses the most popular QI, without considering any other
+// values, a ViewEdge with FEdges having QIs of 0, 21, 22, 23, 24 and 25
+// will end up having a total QI of 0, even though most of the FEdges are
+// heavily occluded.  computeCumulativeVisibility will treat this case as
+// a QI of 22 because 3 out of 6 occluders have QI <= 22.
+
+template <typename G, typename I>
+static void computeCumulativeVisibility(ViewMap *ioViewMap, G& grid, real epsilon)
+{
+  vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
+
+  FEdge * fe, *festart;
+  int nSamples = 0;
+  vector<WFace*> wFaces;
+  WFace *wFace = 0;
+  unsigned tmpQI = 0;
+  unsigned qiClasses[256];
+  unsigned maxIndex, maxCard;
+  unsigned qiMajority;
+  for(vector<ViewEdge*>::iterator ve=vedges.begin(), veend=vedges.end(); ve!=veend; ve++) {
+#if logging > 0
+cout << "Processing ViewEdge " << (*ve)->getId() << endl;
+#endif
+    // Find an edge to test
+    if ( ! (*ve)->isInImage() ) {
+      // This view edge has been proscenium culled
+      (*ve)->setQI(255);
+      (*ve)->setaShape(0);
+#if logging > 0
+cout << "\tCulled." << endl;
+#endif
+      continue;
+    }
+
+    // Test edge
+    festart = (*ve)->fedgeA();
+    fe = (*ve)->fedgeA();
+    qiMajority = 0;
+    do {
+      if ( fe != NULL && fe->isInImage() ) {
+        qiMajority++;
+      }
+      fe = fe->nextEdge();
+    } while (fe && fe != festart);
+
+    if ( qiMajority == 0 ) {
+      // There are no occludable FEdges on this ViewEdge
+      // This should be impossible.
+      cout << "View Edge in viewport without occludable FEdges: " << (*ve)->getId() << endl;
+      // We can recover from this error:
+      // Treat this edge as fully visible with no occludee
+      (*ve)->setQI(0);
+      (*ve)->setaShape(0);
+      continue;
+    } else {
+      ++qiMajority;
+      qiMajority >>= 1;
+    }
+#if logging > 0
+cout << "\tqiMajority: " << qiMajority << endl;
+#endif
+
+    tmpQI = 0;
+    maxIndex = 0;
+    maxCard = 0;
+    nSamples = 0;
+    memset(qiClasses, 0, 256 * sizeof(*qiClasses));
+	set<ViewShape*> foundOccluders;
+
+    fe = (*ve)->fedgeA();
+    do
+    {
+      if ( fe == NULL || ! fe->isInImage() ) {
+        fe = fe->nextEdge();
+        continue;
+      }
+      if((maxCard < qiMajority)) {
+        tmpQI = computeVisibility<G, I>(ioViewMap, fe, grid, epsilon, *ve, &wFace, &foundOccluders); //ARB: change &wFace to wFace and use reference in called function
+#if logging > 0
+cout << "\tFEdge: visibility " << tmpQI << endl;
+#endif
+	
+        //ARB: This is an error condition, not an alert condition.
+        // Some sort of recovery or abort is necessary.
+        if(tmpQI >= 256) {
+          cerr << "Warning: too many occluding levels" << endl;
+          //ARB: Wild guess: instead of aborting or corrupting memory, treat as tmpQI == 255
+          tmpQI = 255;
+        }
+
+        if (++qiClasses[tmpQI] > maxCard) {
+          maxCard = qiClasses[tmpQI];
+          maxIndex = tmpQI;
+        }
+      } else {
+        //ARB: FindOccludee is redundant if ComputeRayCastingVisibility has been called
+        findOccludee<G, I>(fe, grid, epsilon, *ve, &wFace); //ARB: change &wFace to wFace and use reference in called function
+#if logging > 0
+cout << "\tFEdge: occludee only (" << (wFace != NULL ? "found" : "not found") << ")" << endl;
+#endif
+      }
+
+      // Store test results
+      if(wFace) { 
+        vector<Vec3r> vertices;
+        for ( int i = 0, numEdges = wFace->numberOfEdges(); i < numEdges; ++i ) {
+          vertices.push_back(Vec3r(wFace->GetVertex(i)->GetVertex()));
+        }
+        Polygon3r poly(vertices, wFace->GetNormal());
+        poly.userdata = (void *) wFace;
+        fe->setaFace(poly);
+        wFaces.push_back(wFace);
+        fe->setOccludeeEmpty(false);
+#if logging > 0
+cout << "\tFound occludee" << endl;
+#endif
+      } else {
+        fe->setOccludeeEmpty(true);
+      }
+
+      ++nSamples;
+      fe = fe->nextEdge();
+    }
+    while((maxCard < qiMajority) && (0!=fe) && (fe!=festart));
+#if logging > 0
+cout << "\tFinished with " << nSamples << " samples, maxCard = " << maxCard << endl;
+#endif
+
+    // ViewEdge
+    // qi --
+    // Find the minimum value that is >= the majority of the QI
+    for ( unsigned count = 0, i = 0; i < 256; ++i ) {
+        count += qiClasses[i];
+        if ( count >= qiMajority ) {
+            (*ve)->setQI(i);
+            break;
+        }
+    }
+    // occluders --
+    // I would rather not have to go through the effort of creating this
+    // this set and then copying out its contents.  Is there a reason why
+    // ViewEdge::_Occluders cannot be converted to a set<>?
+    for(set<ViewShape*>::iterator o=foundOccluders.begin(), oend=foundOccluders.end(); o!=oend; ++o) {
+      (*ve)->AddOccluder((*o));
+    }
+#if logging > 0
+cout << "\tConclusion: QI = " << maxIndex << ", " << (*ve)->occluders_size() << " occluders." << endl;
+#endif
+    // occludee --
+    if(!wFaces.empty())
+    {
+      if(wFaces.size() <= (float)nSamples/2.f)
+      {
+        (*ve)->setaShape(0);
+      }
+      else
+      {
+        ViewShape *vshape = ioViewMap->viewShape((*wFaces.begin())->GetVertex(0)->shape()->GetId());
+        (*ve)->setaShape(vshape);
+      }
+    }
+    
+    wFaces.clear();
+  }    
+}
+
+template <typename G, typename I>
+static void computeDetailedVisibility(ViewMap *ioViewMap, G& grid, real epsilon)
+{
+  vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
+
+  FEdge * fe, *festart;
+  int nSamples = 0;
+  vector<WFace*> wFaces;
+  WFace *wFace = 0;
+  unsigned tmpQI = 0;
+  unsigned qiClasses[256];
+  unsigned maxIndex, maxCard;
+  unsigned qiMajority;
+  for(vector<ViewEdge*>::iterator ve=vedges.begin(), veend=vedges.end(); ve!=veend; ve++) {
+#if logging > 0
+cout << "Processing ViewEdge " << (*ve)->getId() << endl;
+#endif
+    // Find an edge to test
+    if ( ! (*ve)->isInImage() ) {
+      // This view edge has been proscenium culled
+      (*ve)->setQI(255);
+      (*ve)->setaShape(0);
+#if logging > 0
+cout << "\tCulled." << endl;
+#endif
+      continue;
+    }
+
+    // Test edge
+    festart = (*ve)->fedgeA();
+    fe = (*ve)->fedgeA();
+    qiMajority = 0;
+    do {
+      if ( fe != NULL && fe->isInImage() ) {
+        qiMajority++;
+      }
+      fe = fe->nextEdge();
+    } while (fe && fe != festart);
+
+    if ( qiMajority == 0 ) {
+      // There are no occludable FEdges on this ViewEdge
+      // This should be impossible.
+      cout << "View Edge in viewport without occludable FEdges: " << (*ve)->getId() << endl;
+      // We can recover from this error:
+      // Treat this edge as fully visible with no occludee
+      (*ve)->setQI(0);
+      (*ve)->setaShape(0);
+      continue;
+    } else {
+      ++qiMajority;
+      qiMajority >>= 1;
+    }
+#if logging > 0
+cout << "\tqiMajority: " << qiMajority << endl;
+#endif
+
+    tmpQI = 0;
+    maxIndex = 0;
+    maxCard = 0;
+    nSamples = 0;
+    memset(qiClasses, 0, 256 * sizeof(*qiClasses));
+	set<ViewShape*> foundOccluders;
+
+    fe = (*ve)->fedgeA();
+    do
+    {
+      if ( fe == NULL || ! fe->isInImage() ) {
+        fe = fe->nextEdge();
+        continue;
+      }
+      if((maxCard < qiMajority)) {
+        tmpQI = computeVisibility<G, I>(ioViewMap, fe, grid, epsilon, *ve, &wFace, &foundOccluders); //ARB: change &wFace to wFace and use reference in called function
+#if logging > 0
+cout << "\tFEdge: visibility " << tmpQI << endl;
+#endif
+	
+        //ARB: This is an error condition, not an alert condition.
+        // Some sort of recovery or abort is necessary.
+        if(tmpQI >= 256) {
+          cerr << "Warning: too many occluding levels" << endl;
+          //ARB: Wild guess: instead of aborting or corrupting memory, treat as tmpQI == 255
+          tmpQI = 255;
+        }
+
+        if (++qiClasses[tmpQI] > maxCard) {
+          maxCard = qiClasses[tmpQI];
+          maxIndex = tmpQI;
+        }
+      } else {
+        //ARB: FindOccludee is redundant if ComputeRayCastingVisibility has been called
+        findOccludee<G, I>(fe, grid, epsilon, *ve, &wFace); //ARB: change &wFace to wFace and use reference in called function
+#if logging > 0
+cout << "\tFEdge: occludee only (" << (wFace != NULL ? "found" : "not found") << ")" << endl;
+#endif
+      }
+
+      // Store test results
+      if(wFace) { 
+        vector<Vec3r> vertices;
+        for ( int i = 0, numEdges = wFace->numberOfEdges(); i < numEdges; ++i ) {
+          vertices.push_back(Vec3r(wFace->GetVertex(i)->GetVertex()));
+        }
+        Polygon3r poly(vertices, wFace->GetNormal());
+        poly.userdata = (void *) wFace;
+        fe->setaFace(poly);
+        wFaces.push_back(wFace);
+        fe->setOccludeeEmpty(false);
+#if logging > 0
+cout << "\tFound occludee" << endl;
+#endif
+      } else {
+        fe->setOccludeeEmpty(true);
+      }
+
+      ++nSamples;
+      fe = fe->nextEdge();
+    }
+    while((maxCard < qiMajority) && (0!=fe) && (fe!=festart));
+#if logging > 0
+cout << "\tFinished with " << nSamples << " samples, maxCard = " << maxCard << endl;
+#endif
+
+    // ViewEdge
+    // qi --
+    (*ve)->setQI(maxIndex);
+    // occluders --
+    // I would rather not have to go through the effort of creating this
+    // this set and then copying out its contents.  Is there a reason why
+    // ViewEdge::_Occluders cannot be converted to a set<>?
+    for(set<ViewShape*>::iterator o=foundOccluders.begin(), oend=foundOccluders.end(); o!=oend; ++o) {
+      (*ve)->AddOccluder((*o));
+    }
+#if logging > 0
+cout << "\tConclusion: QI = " << maxIndex << ", " << (*ve)->occluders_size() << " occluders." << endl;
+#endif
+    // occludee --
+    if(!wFaces.empty())
+    {
+      if(wFaces.size() <= (float)nSamples/2.f)
+      {
+        (*ve)->setaShape(0);
+      }
+      else
+      {
+        ViewShape *vshape = ioViewMap->viewShape((*wFaces.begin())->GetVertex(0)->shape()->GetId());
+        (*ve)->setaShape(vshape);
+      }
+    }
+    
+    wFaces.clear();
+  }    
+}
+
+template <typename G, typename I>
+static void computeFastVisibility(ViewMap *ioViewMap, G& grid, real epsilon)
+{
+  vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
+
+  FEdge * fe, *festart;
+  unsigned nSamples = 0;
+  vector<WFace*> wFaces;
+  WFace *wFace = 0;
+  unsigned tmpQI = 0;
+  unsigned qiClasses[256];
+  unsigned maxIndex, maxCard;
+  unsigned qiMajority;
+  bool even_test;
+  for(vector<ViewEdge*>::iterator ve=vedges.begin(), veend=vedges.end(); ve!=veend; ve++) {
+    // Find an edge to test
+    if ( ! (*ve)->isInImage() ) {
+      // This view edge has been proscenium culled
+      (*ve)->setQI(255);
+      (*ve)->setaShape(0);
+      continue;
+    }
+
+    // Test edge
+    festart = (*ve)->fedgeA();
+    fe = (*ve)->fedgeA();
+
+    even_test = true;
+    qiMajority = 0;
+    do {
+      if ( even_test && fe != NULL && fe->isInImage() ) {
+        qiMajority++;
+        even_test = ! even_test;
+      }
+      fe = fe->nextEdge();
+    } while (fe && fe != festart);
+
+    if (qiMajority == 0 ) {
+      // There are no occludable FEdges on this ViewEdge
+      // This should be impossible.
+      cout << "View Edge in viewport without occludable FEdges: " << (*ve)->getId() << endl;
+      // We can recover from this error:
+      // Treat this edge as fully visible with no occludee
+      (*ve)->setQI(0);
+      (*ve)->setaShape(0);
+      continue;
+    } else {
+      ++qiMajority;
+      qiMajority >>= 1;
+    }
+
+    even_test = true;
+    maxIndex = 0;
+    maxCard = 0;
+    nSamples = 0;
+    memset(qiClasses, 0, 256 * sizeof(*qiClasses));
+	set<ViewShape*> foundOccluders;
+
+    fe = (*ve)->fedgeA();
+    do
+    {
+      if ( fe == NULL || ! fe->isInImage() ) {
+        fe = fe->nextEdge();
+        continue;
+      }
+      if (even_test)
+      {
+        if((maxCard < qiMajority)) {
+          tmpQI = computeVisibility<G, I>(ioViewMap, fe, grid, epsilon, *ve, &wFace, &foundOccluders); //ARB: change &wFace to wFace and use reference in called function
+
+          //ARB: This is an error condition, not an alert condition.
+          // Some sort of recovery or abort is necessary.
+          if(tmpQI >= 256) {
+            cerr << "Warning: too many occluding levels" << endl;
+            //ARB: Wild guess: instead of aborting or corrupting memory, treat as tmpQI == 255
+            tmpQI = 255;
+          }
+
+          if (++qiClasses[tmpQI] > maxCard) {
+            maxCard = qiClasses[tmpQI];
+            maxIndex = tmpQI;
+          }
+        }  else {
+          //ARB: FindOccludee is redundant if ComputeRayCastingVisibility has been called
+          findOccludee<G, I>(fe, grid, epsilon, *ve, &wFace); //ARB: change &wFace to wFace and use reference in called function
+        }
+
+        if(wFace)
+        { 
+          vector<Vec3r> vertices;
+          for ( int i = 0, numEdges = wFace->numberOfEdges(); i < numEdges; ++i ) {
+            vertices.push_back(Vec3r(wFace->GetVertex(i)->GetVertex()));
+          }
+          Polygon3r poly(vertices, wFace->GetNormal());
+          poly.userdata = (void *) wFace;
+          fe->setaFace(poly);
+          wFaces.push_back(wFace);
+        } 
+        ++nSamples;
+      }
+
+      even_test = ! even_test;
+      fe = fe->nextEdge();
+    } while ((maxCard < qiMajority) && (0!=fe) && (fe!=festart));
+
+	// qi --
+    (*ve)->setQI(maxIndex);
+
+    // occluders --
+    for(set<ViewShape*>::iterator o=foundOccluders.begin(), oend=foundOccluders.end(); o!=oend; ++o) {
+      (*ve)->AddOccluder((*o));
+    }
+
+    // occludee --
+    if(!wFaces.empty())
+    {
+      if(wFaces.size() < nSamples / 2)
+      {
+        (*ve)->setaShape(0);
+      }
+      else
+      {
+        ViewShape *vshape = ioViewMap->viewShape((*wFaces.begin())->GetVertex(0)->shape()->GetId());
+        (*ve)->setaShape(vshape);
+      }
+    }
+    
+    wFaces.clear();
+  }
+}
+
+template <typename G, typename I>
+static void computeVeryFastVisibility(ViewMap *ioViewMap, G& grid, real epsilon)
+{
+	vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
+
+	FEdge* fe;
+	unsigned qi = 0;
+	WFace* wFace = 0;
+
+	for(vector<ViewEdge*>::iterator ve=vedges.begin(), veend=vedges.end(); ve!=veend; ve++)
+	{
+		// Find an edge to test
+		if ( ! (*ve)->isInImage() ) {
+			// This view edge has been proscenium culled
+			(*ve)->setQI(255);
+			(*ve)->setaShape(0);
+			continue;
+		}
+		fe = (*ve)->fedgeA();
+		// Find a FEdge inside the occluder proscenium to test for visibility
+		FEdge* festart = fe;
+		while ( fe != NULL && ! fe->isInImage() ) {
+			fe = fe->nextEdge();
+			if ( fe == festart ) {
+				break;
+			}
+		}
+
+		// Test edge
+		if ( fe == NULL || ! fe->isInImage() ) {
+			// There are no occludable FEdges on this ViewEdge
+			// This should be impossible.
+			cout << "View Edge in viewport without occludable FEdges: " << (*ve)->getId() << endl;
+			// We can recover from this error:
+			// Treat this edge as fully visible with no occludee
+			qi = 0;
+			wFace = NULL;
+		} else {
+			qi = computeVisibility<G, I>(ioViewMap, fe, grid, epsilon, *ve, &wFace, NULL);
+		}
+
+		// Store test results
+		if(wFace)
+		{ 
+			vector<Vec3r> vertices;
+			for ( int i = 0, numEdges = wFace->numberOfEdges(); i < numEdges; ++i ) {
+				vertices.push_back(Vec3r(wFace->GetVertex(i)->GetVertex()));
+			}
+			Polygon3r poly(vertices, wFace->GetNormal());
+			poly.userdata = (void *) wFace;
+			fe->setaFace(poly);  // This works because setaFace *copies* the polygon
+			ViewShape *vshape = ioViewMap->viewShape(wFace->GetVertex(0)->shape()->GetId());
+			(*ve)->setaShape(vshape);
+		} 
+		else
+		{
+			(*ve)->setaShape(0);
+		}
+		(*ve)->setQI(qi);
+	}  
+
+}
+
+void ViewMapBuilder::BuildGrid(WingedEdge& we, const BBox<Vec3r>& bbox, unsigned int sceneNumFaces) {
+  _Grid->clear();
+  Vec3r size;
+  for(unsigned int i=0; i<3; i++)
+    {
+      size[i] = fabs(bbox.getMax()[i] - bbox.getMin()[i]);
+      size[i] += size[i]/10.0; // let make the grid 1/10 bigger to avoid numerical errors while computing triangles/cells intersections
+      if(size[i]==0){
+          cout << "Warning: the bbox size is 0 in dimension "<<i<<endl;
+      }
+    }
+  _Grid->configure(Vec3r(bbox.getMin() - size / 20.0), size, sceneNumFaces);
+
+  // Fill in the grid:
+  WFillGrid fillGridRenderer(_Grid, &we);
+  fillGridRenderer.fillGrid();
+  
+  // DEBUG
+  _Grid->displayDebug();
+}
+
+ViewMap* ViewMapBuilder::BuildViewMap(WingedEdge& we, visibility_algo iAlgo, real epsilon, 
+      const BBox<Vec3r>& bbox, unsigned int sceneNumFaces) {
   _ViewMap = new ViewMap;
   _currentId = 1;
   _currentFId = 0;
   _currentSVertexId = 0;
-  
+
   // Builds initial view edges
   computeInitialViewEdges(we);
   
@@ -40,11 +922,191 @@ ViewMap* ViewMapBuilder::BuildViewMap(WingedEdge& we, visibility_algo iAlgo, rea
   ComputeIntersections(_ViewMap, sweep_line, epsilon);
   
   // Compute visibility
-  ComputeEdgesVisibility(_ViewMap, iAlgo, _Grid, epsilon);
-  
+  ComputeEdgesVisibility(_ViewMap, we, bbox, sceneNumFaces, iAlgo, epsilon);
+
   return _ViewMap;
 }
 
+static inline real distance2D(const Vec3r & point, const real origin[2]) {
+	return ::hypot((point[0] - origin[0]), (point[1] - origin[1]));
+}
+
+static inline bool crossesProscenium(real proscenium[4], FEdge *fe) {
+	Vec2r min(proscenium[0], proscenium[2]);
+	Vec2r max(proscenium[1], proscenium[3]);
+	Vec2r A(fe->vertexA()->getProjectedX(), fe->vertexA()->getProjectedY());
+	Vec2r B(fe->vertexB()->getProjectedX(), fe->vertexB()->getProjectedY());
+
+	return GeomUtils::intersect2dSeg2dArea (min, max, A, B);
+}
+
+static inline bool insideProscenium(real proscenium[4], const Vec3r& point) {
+	return ! ( point[0] < proscenium[0] || point[0] > proscenium[1] || point[1] < proscenium[2] || point[1] > proscenium[3] );
+}
+
+void ViewMapBuilder::CullViewEdges(ViewMap *ioViewMap, real viewProscenium[4], real occluderProscenium[4], bool extensiveFEdgeSearch) {
+	// Cull view edges by marking them as non-displayable.
+	// This avoids the complications of trying to delete
+	// edges from the ViewMap.
+
+	// Non-displayable view edges will be skipped over during
+	// visibility calculation.
+
+	// View edges will be culled according to their position
+	// w.r.t. the viewport proscenium (viewport + 5% border,
+	// or some such).
+
+	// Get proscenium boundary for culling
+	GridHelpers::getDefaultViewProscenium(viewProscenium);
+	real prosceniumOrigin[2];
+	prosceniumOrigin[0] = (viewProscenium[1] - viewProscenium[0]) / 2.0;
+	prosceniumOrigin[1] = (viewProscenium[3] - viewProscenium[2]) / 2.0;
+	cout << "Proscenium culling:" << endl;
+	cout << "Proscenium: [" << viewProscenium[0] << ", " << viewProscenium[1] << ", " << viewProscenium[2] << ", " << viewProscenium[3] << "]"<< endl;
+	cout << "Origin: [" << prosceniumOrigin[0] << ", " << prosceniumOrigin[1] << "]"<< endl;
+	
+	// A separate occluder proscenium will also be maintained,
+	// starting out the same as the viewport proscenium, and 
+	// expanding as necessary so that it encompasses the center 
+	// point of at least one feature edge in each retained view
+	// edge.
+	// The occluder proscenium will be used later to cull occluding
+	// triangles before they are inserted into the Grid.
+	// The occluder proscenium starts out the same size as the view
+	// proscenium
+	GridHelpers::getDefaultViewProscenium(occluderProscenium);
+
+	// N.B. Freestyle is inconsistent in its use of ViewMap::viewedges_container
+	// and vector<ViewEdge*>::iterator.  Probably all occurences of vector<ViewEdge*>::iterator
+	// should be replaced ViewMap::viewedges_container throughout the code.
+	// For each view edge
+	ViewMap::viewedges_container::iterator ve, veend;
+
+	for(ve=ioViewMap->ViewEdges().begin(), veend=ioViewMap->ViewEdges().end(); ve!=veend; ve++) {
+		// Overview:
+		//    Search for a visible feature edge
+		//    If none: mark view edge as non-displayable
+		//    Otherwise:
+		//        Find a feature edge with center point inside occluder proscenium.
+		//        If none exists, find the feature edge with center point
+		//            closest to viewport origin.
+		//            Expand occluder proscenium to enclose center point.
+
+		// For each feature edge, while bestOccluderTarget not found and view edge not visibile
+		bool bestOccluderTargetFound = false;
+		FEdge *bestOccluderTarget = NULL;
+		real bestOccluderDistance = 0.0;
+		FEdge *festart = (*ve)->fedgeA();
+		FEdge *fe = festart;
+		// All ViewEdges start culled
+		(*ve)->setIsInImage(false);
+
+		// For simple visibility calculation: mark a feature edge
+		// that is known to have a center point inside the occluder proscenium.  
+		// Cull all other feature edges.
+		do {
+			// All FEdges start culled
+			fe->setIsInImage(false);
+
+			// Look for the visible edge that can most easily be included
+			// in the occluder proscenium.
+			if ( ! bestOccluderTargetFound ) {
+				// If center point is inside occluder proscenium, 
+				if ( insideProscenium(occluderProscenium, fe->center2d()) ) {
+					// Use this feature edge for visibility deterimination
+					fe->setIsInImage(true);
+					// Mark bestOccluderTarget as found
+					bestOccluderTargetFound = true;
+					bestOccluderTarget = fe;
+				} else {
+					real d = distance2D(fe->center2d(), prosceniumOrigin);
+					// If center point is closer to viewport origin than current target
+					if ( bestOccluderTarget == NULL || d < bestOccluderDistance ) {
+						// Then store as bestOccluderTarget
+						bestOccluderDistance = d;
+						bestOccluderTarget = fe;
+					}
+				}
+			}
+
+			// If feature edge crosses the view proscenium
+			if ( ! (*ve)->isInImage() && crossesProscenium(viewProscenium, fe) ) {
+				// Then the view edge will be included in the image
+				(*ve)->setIsInImage(true);
+			}
+			fe = fe->nextEdge();
+		} while ( fe != NULL && fe != festart && ! ( bestOccluderTargetFound && (*ve)->isInImage() ) );
+
+		// Either we have run out of FEdges, or we already have the one edge we need to determine visibility
+		// Cull all remaining edges.
+		while ( fe != NULL && fe != festart ) {
+			fe->setIsInImage(false);
+			fe = fe->nextEdge();
+		}
+
+		// If bestOccluderTarget was not found inside the occluder proscenium, 
+		// we need to expand the occluder proscenium to include it.
+		if ( (*ve)->isInImage() && bestOccluderTarget != NULL && ! bestOccluderTargetFound ) {
+			// Expand occluder proscenium to enclose bestOccluderTarget
+			Vec3r point = bestOccluderTarget->center2d();
+			if ( point[0] < occluderProscenium[0] ) {
+				occluderProscenium[0] = point[0];
+			} else if ( point[0] > occluderProscenium[1] ) {
+				occluderProscenium[1] = point[0];
+			}
+			if ( point[1] < occluderProscenium[2] ) {
+				occluderProscenium[2] = point[1];
+			} else if ( point[1] > occluderProscenium[3] ) {
+				occluderProscenium[3] = point[1];
+			}
+			// Use bestOccluderTarget for visibility determination
+			bestOccluderTarget->setIsInImage(true);
+		}
+	}
+
+	// We are done calculating the occluder proscenium.
+	// Expand the occluder proscenium by an epsilon to avoid rounding errors.
+	const real epsilon = 1.0e-6;
+	occluderProscenium[0] -= epsilon;
+	occluderProscenium[1] += epsilon;
+	occluderProscenium[2] -= epsilon;
+	occluderProscenium[3] += epsilon;
+
+	// For "Normal" or "Fast" style visibility computation only:
+
+	// For more detailed visibility calculation, make a second pass through
+	// the view map, marking all feature edges with center points inside
+	// the final occluder proscenium.  All of these feature edges can be
+	// considered during visibility calculation.
+
+	// So far we have only found one FEdge per ViewEdge.  The "Normal" and
+	// "Fast" styles of visibility computation want to consider many
+	// FEdges for each ViewEdge.
+	// Here we re-scan the view map to find any usable FEdges that we
+	// skipped on the first pass, or that have become usable because the
+	// occluder proscenium has been expanded since the edge was visited
+	// on the first pass.
+	if ( extensiveFEdgeSearch ) {
+		// For each view edge,
+		for(ve=ioViewMap->ViewEdges().begin(), veend=ioViewMap->ViewEdges().end(); ve!=veend; ve++) {
+			if ( ! (*ve)->isInImage() ) {
+				continue;
+			}
+			// For each feature edge,
+			FEdge *festart = (*ve)->fedgeA();
+			FEdge *fe = festart;
+			do {
+				// If not (already) visible and center point inside occluder proscenium, 
+				if ( ! fe->isInImage() && insideProscenium(occluderProscenium, fe->center2d()) ) {
+					// Use the feature edge for visibility determination
+					fe->setIsInImage(true);
+				}
+				fe = fe->nextEdge();
+			} while ( fe != NULL && fe != festart );
+		}
+	}
+}
+
 void ViewMapBuilder::computeInitialViewEdges(WingedEdge& we)
 {
   vector<WShape*> wshapes = we.getWShapes();
@@ -151,26 +1213,130 @@ void ViewMapBuilder::computeCusps(ViewMap *ioViewMap){
     vedges.push_back(*ve);
   }
 }
-void ViewMapBuilder::ComputeEdgesVisibility(ViewMap *ioViewMap, visibility_algo iAlgo,  Grid *iGrid, real epsilon)
+
+void ViewMapBuilder::ComputeCumulativeVisibility(ViewMap *ioViewMap, 
+      WingedEdge& we, const BBox<Vec3r>& bbox, real epsilon, bool cull, GridDensityProviderFactory& factory)
 {
-  if((iAlgo == ray_casting ||
-      iAlgo == ray_casting_fast ||
-      iAlgo == ray_casting_very_fast) && (NULL == iGrid))
-  {
-    cerr << "Error: can't cast ray, no grid defined" << endl;
-    return;
-  }
+	auto_ptr<GridHelpers::Transform> transform;
+	auto_ptr<OccluderSource> source;
 
+	if ( _orthographicProjection ) {
+		transform.reset(new BoxGrid::Transform);
+	} else {
+		transform.reset(new SphericalGrid::Transform);
+	}
+
+	if ( cull ) {
+		source.reset(new CulledOccluderSource(*transform, we, *ioViewMap, true));
+	} else {
+		source.reset(new OccluderSource(*transform, we));
+	}
+
+	auto_ptr<GridDensityProvider> density(factory.newGridDensityProvider(*source, bbox, *transform));
+
+	if ( _orthographicProjection ) {
+		BoxGrid grid(*source, *density, ioViewMap, _viewpoint, _EnableQI);
+		computeCumulativeVisibility<BoxGrid, BoxGrid::Iterator>(ioViewMap, grid, epsilon);
+	} else {
+		SphericalGrid grid(*source, *density, ioViewMap, _viewpoint, _EnableQI);
+		computeCumulativeVisibility<SphericalGrid, SphericalGrid::Iterator>(ioViewMap, grid, epsilon);
+	}
+}
+
+void ViewMapBuilder::ComputeDetailedVisibility(ViewMap *ioViewMap, 
+      WingedEdge& we, const BBox<Vec3r>& bbox, real epsilon, bool cull, GridDensityProviderFactory& factory)
+{
+	auto_ptr<GridHelpers::Transform> transform;
+	auto_ptr<OccluderSource> source;
+
+	if ( _orthographicProjection ) {
+		transform.reset(new BoxGrid::Transform);
+	} else {
+		transform.reset(new SphericalGrid::Transform);
+	}
+
+	if ( cull ) {
+		source.reset(new CulledOccluderSource(*transform, we, *ioViewMap, true));
+	} else {
+		source.reset(new OccluderSource(*transform, we));
+	}
+
+	auto_ptr<GridDensityProvider> density(factory.newGridDensityProvider(*source, bbox, *transform));
+
+	if ( _orthographicProjection ) {
+		BoxGrid grid(*source, *density, ioViewMap, _viewpoint, _EnableQI);
+		computeDetailedVisibility<BoxGrid, BoxGrid::Iterator>(ioViewMap, grid, epsilon);
+	} else {
+		SphericalGrid grid(*source, *density, ioViewMap, _viewpoint, _EnableQI);
+		computeDetailedVisibility<SphericalGrid, SphericalGrid::Iterator>(ioViewMap, grid, epsilon);
+	}
+}
+
+void ViewMapBuilder::ComputeEdgesVisibility(ViewMap *ioViewMap, 
+      WingedEdge& we, const BBox<Vec3r>& bbox, unsigned int sceneNumFaces, visibility_algo iAlgo, real epsilon)
+{
   switch(iAlgo)
   {
   case ray_casting:
-    ComputeRayCastingVisibility(ioViewMap, iGrid, epsilon);
+    cout << "Using ordinary ray casting" << endl;
+    BuildGrid(we, bbox, sceneNumFaces);
+    ComputeRayCastingVisibility(ioViewMap, epsilon);
     break;
   case ray_casting_fast:
-    ComputeFastRayCastingVisibility(ioViewMap, iGrid, epsilon);
+    cout << "Using fast ray casting" << endl;
+    BuildGrid(we, bbox, sceneNumFaces);
+    ComputeFastRayCastingVisibility(ioViewMap, epsilon);
     break;
   case ray_casting_very_fast:
-    ComputeVeryFastRayCastingVisibility(ioViewMap, iGrid, epsilon);
+    cout << "Using very fast ray casting" << endl;
+    BuildGrid(we, bbox, sceneNumFaces);
+    ComputeVeryFastRayCastingVisibility(ioViewMap, epsilon);
+    break;
+  case ray_casting_culled_adaptive_traditional:
+	cout << "Using culled adaptive grid with heuristic density and traditional QI calculation" << endl;
+	try {
+		HeuristicGridDensityProviderFactory factory(0.5f, sceneNumFaces);
+		ComputeDetailedVisibility(ioViewMap, we, bbox, epsilon, true, factory);
+	} catch (...) {
+		// Last resort catch to make sure RAII semantics hold for OptimizedGrid
+		// Can be replaced with try...catch block around main() if the program
+		// as a whole is converted to RAII
+
+		// This is the little-mentioned caveat of RAII: RAII does not work unless
+		// destructors are always called, but destructors are only called if all
+		// exceptions are caught (or std::terminate() is replaced).
+
+		// We don't actually handle the exception here, so re-throw it
+		// now that our destructors have had a chance to run.
+		throw;
+	}
+    break;
+  case ray_casting_adaptive_traditional:
+	cout << "Using unculled adaptive grid with heuristic density and traditional QI calculation" << endl;
+	try {
+		HeuristicGridDensityProviderFactory factory(0.5f, sceneNumFaces);
+		ComputeDetailedVisibility(ioViewMap, we, bbox, epsilon, false, factory);
+	} catch (...) {
+		throw;
+	}
+    break;
+  case ray_casting_culled_adaptive_cumulative:
+	cout << "Using culled adaptive grid with heuristic density and cumulative QI calculation" << endl;
+	try {
+		HeuristicGridDensityProviderFactory factory(0.5f, sceneNumFaces);
+		ComputeCumulativeVisibility(ioViewMap, we, bbox, epsilon, true, factory);
+	} catch (...) {
+		throw;
+	}
+    break;
+  case ray_casting_adaptive_cumulative:
+	cout << "Using unculled adaptive grid with heuristic density and cumulative QI calculation" << endl;
+	try {
+		HeuristicGridDensityProviderFactory factory(0.5f, sceneNumFaces);
+		ComputeCumulativeVisibility(ioViewMap, we, bbox, epsilon, false, factory);
+	} catch (...) {
+		throw;
+	}
     break;
   default:
     break;
@@ -180,7 +1346,7 @@ void ViewMapBuilder::ComputeEdgesVisibility(ViewMap *ioViewMap, visibility_algo
 static const unsigned gProgressBarMaxSteps = 10;
 static const unsigned gProgressBarMinSize = 2000;
 
-void ViewMapBuilder::ComputeRayCastingVisibility(ViewMap *ioViewMap, Grid* iGrid, real epsilon)
+void ViewMapBuilder::ComputeRayCastingVisibility(ViewMap *ioViewMap, real epsilon)
 {
   vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
   bool progressBarDisplay = false;
@@ -212,6 +1378,9 @@ void ViewMapBuilder::ComputeRayCastingVisibility(ViewMap *ioViewMap, Grid* iGrid
   ve!=veend;
   ve++)
   {
+#if logging > 0
+cout << "Processing ViewEdge " << (*ve)->getId() << endl;
+#endif
     festart = (*ve)->fedgeA();
     fe = (*ve)->fedgeA();
     qiMajority = 1;
@@ -220,6 +1389,9 @@ void ViewMapBuilder::ComputeRayCastingVisibility(ViewMap *ioViewMap, Grid* iGrid
        fe = fe->nextEdge();
     } while (fe && fe != festart);
     qiMajority >>= 1;
+#if logging > 0
+cout << "\tqiMajority: " << qiMajority << endl;
+#endif
 
     tmpQI = 0;
     maxIndex = 0;
@@ -231,30 +1403,55 @@ void ViewMapBuilder::ComputeRayCastingVisibility(ViewMap *ioViewMap, Grid* iGrid
     do
     {
       if((maxCard < qiMajority)) {
-	tmpQI = ComputeRayCastingVisibility(fe, iGrid, epsilon, occluders, &aFace, timestamp++);
+	tmpQI = ComputeRayCastingVisibility(fe, _Grid, epsilon, occluders, &aFace, timestamp++);
 	
-	if(tmpQI >= 256)
+#if logging > 0
+cout << "\tFEdge: visibility " << tmpQI << endl;
+#endif
+	//ARB: This is an error condition, not an alert condition.
+	// Some sort of recovery or abort is necessary.
+	if(tmpQI >= 256) {
 	  cerr << "Warning: too many occluding levels" << endl;
+	  //ARB: Wild guess: instead of aborting or corrupting memory, treat as tmpQI == 255
+	  tmpQI = 255;
+	}
 
 	if (++qiClasses[tmpQI] > maxCard) {
 	  maxCard = qiClasses[tmpQI];
 	  maxIndex = tmpQI;
 	}
+      } else {
+	//ARB: FindOccludee is redundant if ComputeRayCastingVisibility has been called
+	FindOccludee(fe, _Grid, epsilon, &aFace, timestamp++);
+#if logging > 0
+cout << "\tFEdge: occludee only (" << (aFace != NULL ? "found" : "not found") << ")" << endl;
+#endif
       }
-      FindOccludee(fe, iGrid, epsilon, &aFace, timestamp++);
 
       if(aFace) { 
 	fe->setaFace(*aFace);
 	aFaces.push_back(aFace);
 	fe->setOccludeeEmpty(false);
+#if logging > 0
+cout << "\tFound occludee" << endl;
+#endif
       }
       else
+      {
+	//ARB: We are arbitrarily using the last observed value for occludee 
+	// (almost always the value observed for the edge before festart).
+	// Is that meaningful?
+	// ...in fact, _occludeeEmpty seems to be unused.
 	fe->setOccludeeEmpty(true);
+      }
 
       ++nSamples;
       fe = fe->nextEdge();
     }
     while((maxCard < qiMajority) && (0!=fe) && (fe!=festart));
+#if logging > 0
+cout << "\tFinished with " << nSamples << " samples, maxCard = " << maxCard << endl;
+#endif
 
     // ViewEdge
     // qi --
@@ -264,6 +1461,9 @@ void ViewMapBuilder::ComputeRayCastingVisibility(ViewMap *ioViewMap, Grid* iGrid
     o!=oend;
     ++o)
       (*ve)->AddOccluder((*o));
+#if logging > 0
+cout << "\tConclusion: QI = " << maxIndex << ", " << (*ve)->occluders_size() << " occluders." << endl;
+#endif
     // occludee --
     if(!aFaces.empty())
     {
@@ -292,7 +1492,7 @@ void ViewMapBuilder::ComputeRayCastingVisibility(ViewMap *ioViewMap, Grid* iGrid
   }    
 }
 
-void ViewMapBuilder::ComputeFastRayCastingVisibility(ViewMap *ioViewMap, Grid* iGrid, real epsilon)
+void ViewMapBuilder::ComputeFastRayCastingVisibility(ViewMap *ioViewMap, real epsilon)
 {
   vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
   bool progressBarDisplay = false;
@@ -350,17 +1550,24 @@ void ViewMapBuilder::ComputeFastRayCastingVisibility(ViewMap *ioViewMap, Grid* i
       if (even_test)
       {
 	if((maxCard < qiMajority)) {
-	  tmpQI = ComputeRayCastingVisibility(fe, iGrid, epsilon, occluders, &aFace, timestamp++);
+	  tmpQI = ComputeRayCastingVisibility(fe, _Grid, epsilon, occluders, &aFace, timestamp++);
 	  
-	  if(tmpQI >= 256)
-	    cerr << "Warning: too many occluding levels" << endl;
+	//ARB: This is an error condition, not an alert condition.
+	// Some sort of recovery or abort is necessary.
+	if(tmpQI >= 256) {
+	  cerr << "Warning: too many occluding levels" << endl;
+	  //ARB: Wild guess: instead of aborting or corrupting memory, treat as tmpQI == 255
+	  tmpQI = 255;
+	}
 
 	  if (++qiClasses[tmpQI] > maxCard) {
 	    maxCard = qiClasses[tmpQI];
 	    maxIndex = tmpQI;
 	  }
-	}
-        FindOccludee(fe, iGrid, epsilon, &aFace, timestamp++);
+      }  else {
+		//ARB: FindOccludee is redundant if ComputeRayCastingVisibility has been called
+	       	FindOccludee(fe, _Grid, epsilon, &aFace, timestamp++);
+      }
 
         if(aFace)
         { 
@@ -419,7 +1626,7 @@ void ViewMapBuilder::ComputeFastRayCastingVisibility(ViewMap *ioViewMap, Grid* i
   }
 }
 
-void ViewMapBuilder::ComputeVeryFastRayCastingVisibility(ViewMap *ioViewMap, Grid* iGrid, real epsilon)
+void ViewMapBuilder::ComputeVeryFastRayCastingVisibility(ViewMap *ioViewMap, real epsilon)
 {
   vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
   bool progressBarDisplay = false;
@@ -449,7 +1656,7 @@ void ViewMapBuilder::ComputeVeryFastRayCastingVisibility(ViewMap *ioViewMap, Gri
     set<ViewShape*> occluders;
 
     fe = (*ve)->fedgeA();
-    qi = ComputeRayCastingVisibility(fe, iGrid, epsilon, occluders, &aFace, timestamp++);
+    qi = ComputeRayCastingVisibility(fe, _Grid, epsilon, occluders, &aFace, timestamp++);
     if(aFace)
     { 
       fe->setaFace(*aFace);
@@ -474,7 +1681,6 @@ void ViewMapBuilder::ComputeVeryFastRayCastingVisibility(ViewMap *ioViewMap, Gri
   }  
 }
 
-
 void ViewMapBuilder::FindOccludee(FEdge *fe, Grid* iGrid, real epsilon, Polygon3r** oaPolygon, unsigned timestamp, 
 				  Vec3r& u, Vec3r& A, Vec3r& origin, Vec3r& edge, vector<WVertex*>& faceVertices)
 {
@@ -685,17 +1891,31 @@ int ViewMapBuilder::ComputeRayCastingVisibility(FEdge *fe, Grid* iGrid, real eps
       // this edge).
       //-----------
       oface = (WFace*)(*p)->userdata;
+#if logging > 1
+cout << "\t\tEvaluating intersection for occluder " << ((*p)->getVertices())[0] << ((*p)->getVertices())[1] << ((*p)->getVertices())[2] << endl << "\t\t\tand ray " << vp << " * " << u << " (center " << center << ")" << endl;
+#endif
       Vec3r v1(((*p)->getVertices())[0]);
       Vec3r normal((*p)->getNormal());
       real d = -(v1 * normal);
       real t, t_u, t_v;
             
+#if logging > 1
+cout << "\t\tp:  " << ((*p)->getVertices())[0] << ((*p)->getVertices())[1] << ((*p)->getVertices())[2] << ", norm: " << (*p)->getNormal() << endl;
+#endif
+
       if(0 != face)
 	{
+#if logging > 1
+cout << "\t\tDetermining face adjacency...";
+#endif
 	  skipFace = false;
               
-	  if(face == oface)
+	  if(face == oface) {
+#if logging > 1
+cout << "  Rejecting occluder for face concurrency." << endl;
+#endif
 	    continue;
+      }
               
               
 	  for(vector<WVertex*>::iterator fv=faceVertices.begin(), fvend=faceVertices.end();
@@ -724,8 +1944,12 @@ int ViewMapBuilder::ComputeRayCastingVisibility(FEdge *fe, Grid* iGrid, real eps
 	      if(skipFace)
 		break;
 	    }
-	  if(skipFace)
+	  if(skipFace) {
+#if logging > 1
+cout << "  Rejecting occluder for face adjacency." << endl;
+#endif
 	    continue;
+      }
 	}
       else
 	{
@@ -733,15 +1957,28 @@ int ViewMapBuilder::ComputeRayCastingVisibility(FEdge *fe, Grid* iGrid, real eps
 	  //-------------------------------------------------------------
 	  //first let us compute the plane equation.
            
-	  if(GeomUtils::COINCIDENT == GeomUtils::intersectRayPlane(origin, edge, normal, d, t, epsilon))
+	  if(GeomUtils::COINCIDENT == GeomUtils::intersectRayPlane(origin, edge, normal, d, t, epsilon)) {
+#if logging > 1
+cout << "\t\tRejecting occluder for target coincidence." << endl;
+#endif
 	    continue;
+      }
 	}
 
       if((*p)->rayIntersect(center, u, t, t_u, t_v))
 	{
+#if logging > 1
+cout << "\t\tRay " << vp << " * " << u << " intersects at time " << t << " (raylength is " << raylength << ")" << endl;
+#endif
+#if logging > 1
+cout << "\t\t(u * normal) == " << (u * normal) << " for normal " << normal << endl;
+#endif
 	  if (fabs(u * normal) > 0.0001)
 	    if ((t>0.0) && (t<raylength))
 	      {
+#if logging > 1
+cout << "\t\tIs occluder" << endl;
+#endif
 		WFace *f = (WFace*)((*p)->userdata);
 		ViewShape *vshape = _ViewMap->viewShape(f->GetVertex(0)->shape()->GetId());
 		oOccluders.insert(vshape);
@@ -1051,13 +2288,13 @@ void ViewMapBuilder::ComputeSweepLineIntersections(ViewMap *ioViewMap, real epsi
   fe++)
     (*fe)->userdata = NULL;
 
-  // delete segments
-  //  if(!segments.empty()){
-  //    for(s=segments.begin(),send=segments.end();
-  //    s!=send;
-  //    s++){
-  //      delete *s;
-  //    }
-  //    segments.clear();
-  //  }
+   // delete segments
+   if(!segments.empty()){
+     for(s=segments.begin(),send=segments.end();
+     s!=send;
+     s++){
+       delete *s;
+     }
+   }
 }
+
diff --git a/source/blender/freestyle/intern/view_map/ViewMapBuilder.h b/source/blender/freestyle/intern/view_map/ViewMapBuilder.h
index 12b1266fa12..2e37a41df3c 100755
--- a/source/blender/freestyle/intern/view_map/ViewMapBuilder.h
+++ b/source/blender/freestyle/intern/view_map/ViewMapBuilder.h
@@ -46,7 +46,8 @@
 # include "../scene_graph/TriangleRep.h"
 # include "../winged_edge/WEdge.h"
 # include "ViewEdgeXBuilder.h"
-
+# include "../system/TimeUtils.h"
+# include "GridDensityProvider.h"
 
 using namespace Geometry;
 
@@ -64,6 +65,7 @@ private:
   bool _EnableQI;
   double _epsilon;
 
+
   // tmp values:
   int _currentId;
   int _currentFId;
@@ -79,7 +81,11 @@ public:
   typedef enum {
     ray_casting,
     ray_casting_fast,
-    ray_casting_very_fast
+    ray_casting_very_fast,
+    ray_casting_culled_adaptive_traditional,
+    ray_casting_adaptive_traditional,
+    ray_casting_culled_adaptive_cumulative,
+    ray_casting_adaptive_cumulative
   } visibility_algo;
 
   inline ViewMapBuilder()
@@ -101,6 +107,10 @@ public:
     }
   }
 
+  /* Build Grid for ray casting */
+  /*! Build non-culled Grid in camera space for ray casting */
+  void BuildGrid(WingedEdge& we, const BBox<Vec3r>& bbox, unsigned int sceneNumFaces);
+
   /*! Compute Shapes from a WingedEdge containing a list of WShapes */
   void computeInitialViewEdges(WingedEdge&);
 
@@ -145,7 +155,9 @@ public:
    *      The root group node containing the WEdge structured scene 
    */
 
-  ViewMap* BuildViewMap(WingedEdge& we, visibility_algo iAlgo = ray_casting, real epsilon=1e-06) ;
+  ViewMap* BuildViewMap(WingedEdge& we, visibility_algo iAlgo, real epsilon, 
+      const BBox<Vec3r>& bbox, unsigned int sceneNumFaces);
+  void CullViewEdges(ViewMap *ioViewMap, real viewProscenium[4], real occluderProscenium[4], bool extensiveFEdgeSearch = true);
   /*! computes the intersection between all 2D
    *  feature edges of the scene.
    *    ioViewMap
@@ -164,7 +176,8 @@ public:
    *    iGrid
    *      For the Ray Casting algorithm. 
    */
-  void ComputeEdgesVisibility(ViewMap *ioViewMap, visibility_algo iAlgo= ray_casting, Grid* iGrid = 0, real epsilon=1e-6);
+  void ComputeEdgesVisibility(ViewMap *ioViewMap, WingedEdge& we, const BBox<Vec3r>& bbox, unsigned int sceneNumFaces, 
+      visibility_algo iAlgo= ray_casting, real epsilon=1e-6);
 
   void setGrid(Grid *iGrid) {_Grid = iGrid;}
   
@@ -192,9 +205,14 @@ protected:
    *      The visibility corresponding to each edge of ioScene is set is this
    *      edge.
    */
-  void ComputeRayCastingVisibility(ViewMap *ioViewMap, Grid *iGrid, real epsilon=1e-6);
-  void ComputeFastRayCastingVisibility(ViewMap *ioViewMap, Grid *iGrid, real epsilon=1e-6);
-  void ComputeVeryFastRayCastingVisibility(ViewMap *ioViewMap, Grid *iGrid, real epsilon=1e-6);
+  void ComputeRayCastingVisibility(ViewMap *ioViewMap, real epsilon=1e-6);
+  void ComputeFastRayCastingVisibility(ViewMap *ioViewMap, real epsilon=1e-6);
+  void ComputeVeryFastRayCastingVisibility(ViewMap *ioViewMap, real epsilon=1e-6);
+
+void ComputeCumulativeVisibility(ViewMap *ioViewMap, WingedEdge& we, 
+	const BBox<Vec3r>& bbox, real epsilon, bool cull, GridDensityProviderFactory& factory);
+void ComputeDetailedVisibility(ViewMap *ioViewMap, WingedEdge& we, 
+	const BBox<Vec3r>& bbox, real epsilon, bool cull, GridDensityProviderFactory& factory);
 
   /*! Compute the visibility for the FEdge fe.
    *  The occluders are added to fe occluders list.