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.

1 files changed, 1277 insertions, 40 deletions

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;
+     }
+   }
 }
+