diff --git a/simgear/scene/tgdb/ShaderGeometry.cxx b/simgear/scene/tgdb/ShaderGeometry.cxx
index d98b5b66..7fdf657f 100644
--- a/simgear/scene/tgdb/ShaderGeometry.cxx
+++ b/simgear/scene/tgdb/ShaderGeometry.cxx
@@ -1,6 +1,32 @@
+/* -*-c++-*-
+ *
+ * Copyright (C) 2008 Stuart Buchanan
+ *
+ * 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., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ *
+ */
+
+#include <osgDB/Registry>
+#include <osgDB/Input>
+#include <osgDB/ParameterOutput>
+
 #include "ShaderGeometry.hxx"
 
 using namespace osg;
+using namespace osgDB;
 
 namespace simgear
 {
@@ -29,4 +55,69 @@ BoundingBox ShaderGeometry::computeBound() const
     return bb;
 }
 
+bool ShaderGeometry_readLocalData(Object& obj, Input& fr)
+{
+    bool iteratorAdvanced = false;
+
+    ShaderGeometry& geom = static_cast<ShaderGeometry&>(obj);
+
+    if ((fr[0].matchWord("geometry"))) {
+        ++fr;
+        iteratorAdvanced = true;
+        Drawable* drawable = fr.readDrawable();
+        if (drawable) {
+            geom._geometry = drawable;
+        }
+    }
+    if ((fr.matchSequence("instances %i"))) {
+        int entry = fr[0].getNoNestedBrackets();
+        int capacity;
+        fr[1].getInt(capacity);
+        geom._trees.reserve(capacity);
+        fr += 3;
+        iteratorAdvanced = true;
+        // skip {
+        while (!fr.eof() && fr[0].getNoNestedBrackets() > entry) {
+            Vec4 v;
+            if (fr[0].getFloat(v.x()) && fr[1].getFloat(v.y())
+                && fr[2].getFloat(v.z()) && fr[3].getFloat(v.w())) {
+                    fr += 4;
+                    geom._trees.push_back(v);
+            } else {
+                ++fr;
+            }
+        }
+    }
+    return iteratorAdvanced;
+}
+
+bool ShaderGeometry_writeLocalData(const Object& obj, Output& fw)
+{
+    const ShaderGeometry& geom = static_cast<const ShaderGeometry&>(obj);
+
+    fw.indent() << "geometry" << std::endl;
+    fw.writeObject(*geom._geometry);
+    fw.indent() << "instances " << geom._trees.size() << std::endl;
+    fw.indent() << "{" << std::endl;
+    fw.moveIn();
+    for (ShaderGeometry::PositionSizeList::const_iterator iter
+             = geom._trees.begin();
+         iter != geom._trees.end();
+         ++iter) {
+        fw.indent() << iter->x() << " " << iter->y() << " " << iter->z() << " "
+                    << iter->w() << std::endl;
+    }
+    fw.moveOut();
+    fw.indent() << "}" << std::endl;
+    return true;
+}
+
+osgDB::RegisterDotOsgWrapperProxy shaderGeometryProxy
+(
+    new ShaderGeometry,
+    "ShaderGeometry",
+    "Object Drawable ShaderGeometry",
+    &ShaderGeometry_readLocalData,
+    &ShaderGeometry_writeLocalData
+    );
 }
diff --git a/simgear/scene/tgdb/ShaderGeometry.hxx b/simgear/scene/tgdb/ShaderGeometry.hxx
index 78eb1d4f..99f2c1e0 100644
--- a/simgear/scene/tgdb/ShaderGeometry.hxx
+++ b/simgear/scene/tgdb/ShaderGeometry.hxx
@@ -47,7 +47,7 @@ class ShaderGeometry : public osg::Drawable
         ShaderGeometry(const ShaderGeometry& ShaderGeometry,const osg::CopyOp& copyop=osg::CopyOp::SHALLOW_COPY):
             osg::Drawable(ShaderGeometry,copyop) {}
 
-        META_Object(osg,ShaderGeometry);
+        META_Object(flightgear, ShaderGeometry);
         
         typedef std::vector<osg::Vec4> PositionSizeList;
         
@@ -55,7 +55,7 @@ class ShaderGeometry : public osg::Drawable
         virtual osg::BoundingBox computeBound() const;
     
         
-        void setGeometry(osg::Geometry* geometry)
+        void setGeometry(osg::Drawable* geometry)
         {
             _geometry = geometry;
         }
@@ -65,7 +65,7 @@ class ShaderGeometry : public osg::Drawable
             _trees.push_back(osg::Vec4(position, scale));
         }
         
-        osg::ref_ptr<osg::Geometry> _geometry;
+        osg::ref_ptr<osg::Drawable> _geometry;
 
         PositionSizeList _trees;
 
diff --git a/simgear/scene/tgdb/TreeBin.cxx b/simgear/scene/tgdb/TreeBin.cxx
index dfad0dc4..a9777444 100644
--- a/simgear/scene/tgdb/TreeBin.cxx
+++ b/simgear/scene/tgdb/TreeBin.cxx
@@ -19,6 +19,10 @@
  *
  */
 
+#include <algorithm>
+#include <string>
+#include <map>
+
 #include <osg/AlphaFunc>
 #include <osg/Billboard>
 #include <osg/BlendFunc>
@@ -36,11 +40,26 @@
 #include <osgDB/FileUtils>
 
 #include <simgear/misc/sg_path.hxx>
+#include <simgear/scene/util/QuadTreeBuilder.hxx>
+#include <simgear/scene/util/RenderConstants.hxx>
+#include <simgear/scene/util/StateAttributeFactory.hxx>
 
 #include "ShaderGeometry.hxx"
 #include "TreeBin.hxx"
 
-#define SG_TREE_QUAD_TREE_SIZE 32
+#define SG_TREE_QUAD_TREE_DEPTH 3
+
+// Comments from Tim Moore:
+// Some work remains for this code. Stuart's enhancement for multiple
+// textures per forest should be integrated. We should try to use one
+// ShaderGeometry for *all* the trees in the scene graph and do the
+// rotation and scale with a MatrixTransform above the trees quad
+// tree. The positions would of course have to be transformed by the
+// inverse of that transform. Also, we should investigate whether it
+// would be better to instantiate trees as polygons in a osg::Geometry
+// object instead of using the ShaderGeometry instancing technique.
+
+using namespace osg;
 
 namespace simgear
 {
@@ -50,24 +69,12 @@ osg::Geometry* createOrthQuads(float w, float h, const osg::Matrix& rotate)
 
     //const osg::Vec3& pos = osg::Vec3(0.0f,0.0f,0.0f),
     // set up the coords
+    // Create front and back polygons so we don't need to screw around
+    // with two-sided lighting in the shader.
     osg::Vec3Array& v = *(new osg::Vec3Array(8));
+    osg::Vec3Array& n = *(new osg::Vec3Array(8));
     osg::Vec2Array& t = *(new osg::Vec2Array(8));
     
-    /*
-    float rotation = 0.0f;
-    float sw = sinf(rotation)*w*0.5f;
-    float cw = cosf(rotation)*w*0.5f;
-
-    v[0].set(pos.x()-sw,pos.y()-cw,pos.z()+0.0f);
-    v[1].set(pos.x()+sw,pos.y()+cw,pos.z()+0.0f);
-    v[2].set(pos.x()+sw,pos.y()+cw,pos.z()+h);
-    v[3].set(pos.x()-sw,pos.y()-cw,pos.z()+h);
-
-    v[4].set(pos.x()-cw,pos.y()+sw,pos.z()+0.0f);
-    v[5].set(pos.x()+cw,pos.y()-sw,pos.z()+0.0f);
-    v[6].set(pos.x()+cw,pos.y()-sw,pos.z()+h);
-    v[7].set(pos.x()-cw,pos.y()+sw,pos.z()+h);
-    */
     float cw = w*0.5f;
 
     v[0].set(0.0f,-cw,0.0f);
@@ -89,24 +96,120 @@ osg::Geometry* createOrthQuads(float w, float h, const osg::Matrix& rotate)
     t[5].set(1.0f,0.0f);
     t[6].set(1.0f,1.0f);
     t[7].set(0.0f,1.0f);
-    
-    for (unsigned int i = 0; i < 8; i++)
-    {
-      v[i] = v[i] * rotate;      
+
+    // For now the normal is normal to the quad. If we want to get
+    // fancier and approximate a cylindrical tree or something, then
+    // we would really want more geometry.
+    std::fill(n.begin(), n.begin() + 4, Vec3f(1.0f, 0.0f, 0.0f));
+    std::fill(n.begin() + 4, n.end(), Vec3f(0.0f, -1.0f, 0.0f));
+    for (unsigned int i = 0; i < 8; i++) {
+        v[i] = v[i] * rotate;
+        // Should be the inverse transpose, but assume that rotate is
+        // orthonormal.
+        n[i] = n[i] * rotate;     
     }
 
     osg::Geometry *geom = new osg::Geometry;
 
-    geom->setVertexArray( &v );
-
-    geom->setTexCoordArray( 0, &t );
-
-    geom->addPrimitiveSet( new osg::DrawArrays(osg::PrimitiveSet::QUADS,0,8) );
+    geom->setVertexArray(&v);
+    geom->setTexCoordArray(0, &t);
+    geom->setNormalArray(&n);
+    geom->setNormalBinding(Geometry::BIND_PER_VERTEX);
+    // No color for now; that's used to pass the position.
+    geom->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::QUADS,0,8));
 
     return geom;
 }
 
-osg::Group* createForest(const TreeBin& forest, const osg::Matrix& transform)
+ static char vertexShaderSource[] = 
+    "varying float fogFactor;\n"
+    "\n"
+    "void main(void)\n"
+    "{\n"
+    "  vec3 position = gl_Vertex.xyz * gl_Color.w + gl_Color.xyz;\n"
+    "  gl_Position   = gl_ModelViewProjectionMatrix * vec4(position,1.0);\n"
+    "  vec3 ecPosition = vec3(gl_ModelViewMatrix * vec4(position, 1.0));\n"
+    "  vec3 N = normalize(gl_NormalMatrix * gl_Normal);\n"
+    "  vec3 diffuse = gl_FrontMaterial.diffuse.rgb * max(0.0, dot(N, gl_LightSource[0].position.xyz));\n"
+    "  vec3 backDiffuse = gl_FrontMaterial.diffuse.rgb * max(0.0, dot(-N, gl_LightSource[0].position.xyz));\n"
+    " vec4 ambientColor = gl_FrontLightModelProduct.sceneColor + gl_LightSource[0].ambient * gl_FrontMaterial.ambient;\n"
+    " gl_FrontColor = ambientColor + gl_LightSource[0].diffuse * vec4(diffuse, 1.0);\n"
+    " gl_BackColor = ambientColor + gl_LightSource[0].diffuse * vec4(backDiffuse, 1.0)\n;"
+    "  gl_TexCoord[0] = gl_MultiTexCoord0;\n"
+    " float fogCoord = abs(ecPosition.z);\n"
+    "  fogFactor = exp( -gl_Fog.density * gl_Fog.density * fogCoord * fogCoord);\n"
+    "  fogFactor = clamp(fogFactor, 0.0, 1.0);\n"
+    "}\n";
+
+static char fragmentShaderSource[] = 
+    "uniform sampler2D baseTexture; \n"
+//        "varying vec3 N;\n"
+//        "varying vec3 v;\n"
+    "varying float fogFactor;\n"
+    "\n"
+    "void main(void) \n"
+    "{ \n"
+    "  vec4 base = texture2D( baseTexture, gl_TexCoord[0].st);\n"
+    
+    "  vec4 finalColor = base * gl_Color;\n"
+    "  gl_FragColor = mix(gl_Fog.color, finalColor, fogFactor );\n"
+    "}\n";
+
+typedef std::map<std::string, osg::ref_ptr<StateSet> > StateSetMap;
+
+static StateSetMap treeTextureMap;
+
+// Helper classes for creating the quad tree
+namespace
+{
+struct MakeTreesLeaf
+{
+    MakeTreesLeaf(float range, Geometry* geometry) :
+        _range(range), _geometry(geometry)
+    {}
+    MakeTreesLeaf(const MakeTreesLeaf& rhs) :
+        _range(rhs._range), _geometry(rhs._geometry) {}
+    LOD* operator() () const
+    {
+        LOD* result = new LOD;
+        Geode* geode = new Geode;
+        ShaderGeometry* sg = new ShaderGeometry;
+        sg->setGeometry(_geometry);
+        geode->addDrawable(sg);
+        result->addChild(geode, 0, _range);
+        return result;
+    }
+    float _range;
+    Geometry* _geometry;
+};
+
+struct AddTreesLeafObject
+{
+    void operator() (LOD* lod, const TreeBin::Tree& tree) const
+    {
+        Geode* geode = static_cast<Geode*>(lod->getChild(0));
+        ShaderGeometry* sg
+            = static_cast<ShaderGeometry*>(geode->getDrawable(0));
+        sg->addTree(tree.position.osg(), tree.height);
+    }
+};
+
+struct GetTreeCoord
+{
+    GetTreeCoord(const Matrix& transform) : _transform(transform) {}
+    GetTreeCoord(const GetTreeCoord& rhs) : _transform(rhs._transform) {}
+    Vec3 operator() (const TreeBin::Tree& tree) const
+    {
+        return tree.position.osg() * _transform;
+    }
+    Matrix _transform;
+};
+
+typedef QuadTreeBuilder<LOD*, TreeBin::Tree, MakeTreesLeaf, AddTreesLeafObject,
+                        GetTreeCoord> ShaderGeometryQuadtree;
+}
+
+osg::Group* createForest(TreeBin& forest, const osg::Matrix& transform)
 {
     // Set up some shared structures. 
     // FIXME: Currently we only take the texture, height and width of the first tree in the forest. In the future
@@ -116,129 +219,66 @@ osg::Group* createForest(const TreeBin& forest, const osg::Matrix& transform)
     osg::Geometry* shared_geometry = createOrthQuads(firstTree.width, 
                                                      firstTree.height, 
                                                      transform);
-    osg::Group* group = new osg::Group;
+    ref_ptr<Group> group;
 
-    osg::Texture2D *tex = new osg::Texture2D;
-    tex->setWrap( osg::Texture2D::WRAP_S, osg::Texture2D::CLAMP );
-    tex->setWrap( osg::Texture2D::WRAP_T, osg::Texture2D::CLAMP );
-    tex->setImage(osgDB::readImageFile(firstTree.texture));
+    osg::StateSet* stateset = 0;
+    StateSetMap::iterator iter = treeTextureMap.find(firstTree.texture);
+    if (iter == treeTextureMap.end()) {
+        osg::Texture2D *tex = new osg::Texture2D;
+        tex->setWrap( osg::Texture2D::WRAP_S, osg::Texture2D::CLAMP );
+        tex->setWrap( osg::Texture2D::WRAP_T, osg::Texture2D::CLAMP );
+        tex->setImage(osgDB::readImageFile(firstTree.texture));
 
-    osg::AlphaFunc* alphaFunc = new osg::AlphaFunc;
-    alphaFunc->setFunction(osg::AlphaFunc::GEQUAL,0.05f);
-
-    osg::StateSet *dstate = new osg::StateSet;
-    dstate->setTextureAttributeAndModes(0, tex, osg::StateAttribute::ON );
-    dstate->setTextureAttribute(0, new osg::TexEnv );
-    dstate->setAttributeAndModes( new osg::BlendFunc, osg::StateAttribute::ON );
-    dstate->setAttributeAndModes( alphaFunc, osg::StateAttribute::ON );
-    dstate->setMode( GL_LIGHTING, osg::StateAttribute::OFF );
-    dstate->setRenderingHint( osg::StateSet::TRANSPARENT_BIN );
+        static ref_ptr<AlphaFunc> alphaFunc;
+        static ref_ptr<Program> program;
+        static ref_ptr<Uniform> baseTextureSampler;
+        static ref_ptr<Material> material;
     
-    osg::StateSet* stateset = new osg::StateSet;
-    stateset->setTextureAttributeAndModes(0, tex, osg::StateAttribute::ON );
-    stateset->setRenderingHint( osg::StateSet::TRANSPARENT_BIN );
-
-    osg::Program* program = new osg::Program;
-    stateset->setAttribute(program);
-    osg::Uniform* baseTextureSampler = new osg::Uniform("baseTexture",0);
-    stateset->addUniform(baseTextureSampler);
-
-    /*
-     * FIXME: Currently, calculating the diffuse term results in a bad 
-     * "flickering" and a tendency of the diffuse term be either 
-     * 0.0 of 1.0. Hence, it has been commented out in the shader below.
-     * I (Stuart) suspect it may be because the light is so distant that
-     * we're seeing floating point representation issues.
-     */
-    char vertexShaderSource[] = 
-//        "varying vec3 N;\n"
-//        "varying vec3 v;\n"
-        "varying vec2 texcoord;\n"
-        "varying float fogFactor;\n"
-        "\n"
-        "void main(void)\n"
-        "{\n"
-//        "  v = vec3(gl_ModelViewMatrix * gl_Vertex);\n"
-//        "  N = normalize(gl_NormalMatrix * gl_Normal);\n"
-        "  texcoord = gl_MultiTexCoord0.st;\n"        
-        "  vec3 position = gl_Vertex.xyz * gl_Color.w + gl_Color.xyz;\n"
-        "  gl_Position   = gl_ModelViewProjectionMatrix * vec4(position,1.0);\n"
-	    "  const float LOG2 = 1.442695;\n"
-	    "  gl_FogFragCoord = gl_Position.z;\n"
-	    "  fogFactor = exp2( -gl_Fog.density * gl_Fog.density * gl_FogFragCoord * gl_FogFragCoord * LOG2 );\n"
-	    "  fogFactor = clamp(fogFactor, 0.0, 1.0);\n"
-        "}\n";
-
-    char fragmentShaderSource[] = 
-        "uniform sampler2D baseTexture; \n"
-//        "varying vec3 N;\n"
-//        "varying vec3 v;\n"
-        "varying vec2 texcoord;\n"
-        "varying float fogFactor;\n"
-        "\n"
-        "void main(void) \n"
-        "{ \n"
-        "  vec4 base = texture2D( baseTexture, texcoord);\n"
-//        "  vec3 L = normalize(gl_LightSource[0].position.xyz);\n"
-//        "  vec4 vDiffuse = gl_FrontLightProduct[0].diffuse * max(dot(N,L), 0.0);\n"
-//        "  vDiffuse = sqrt(clamp(vDiffuse, 0.0, 1.0));\n"
-//        "  vec4 vAmbient = gl_FrontLightProduct[0].ambient;\n"
-//        "  vec4 finalColor = base * (vAmbient + vDiffuse);\n"
-          "  vec4 finalColor = base * gl_FrontLightProduct[0].diffuse;\n"
-          "  gl_FragColor = mix(gl_Fog.color, finalColor, fogFactor );\n"
-        "}\n";
-
-    osg::Shader* vertex_shader = new osg::Shader(osg::Shader::VERTEX, vertexShaderSource);
-    program->addShader(vertex_shader);
-    
-    osg::Shader* fragment_shader = new osg::Shader(osg::Shader::FRAGMENT, fragmentShaderSource);
-    program->addShader(fragment_shader);
-   
+        stateset = new osg::StateSet;
+        stateset->setTextureAttributeAndModes(0, tex, osg::StateAttribute::ON );
+        stateset->setRenderBinDetails(RANDOM_OBJECTS_BIN, "DepthSortedBin");
+        if (!program.valid()) {
+            alphaFunc = new AlphaFunc;
+            alphaFunc->setFunction(AlphaFunc::GEQUAL,0.33f);
+            program  = new Program;
+            baseTextureSampler = new osg::Uniform("baseTexture", 0);
+            Shader* vertex_shader = new Shader(Shader::VERTEX, vertexShaderSource);
+            program->addShader(vertex_shader);
+            Shader* fragment_shader = new Shader(Shader::FRAGMENT,
+                                                 fragmentShaderSource);
+            program->addShader(fragment_shader);
+            material = new Material;
+            // Don´t track vertex color
+            material->setColorMode(Material::OFF);
+            material->setAmbient(Material::FRONT_AND_BACK,
+                                 Vec4(.6f, .6f, .6f, 1.0f));
+            material->setDiffuse(Material::FRONT_AND_BACK,
+                                 Vec4(.4f, .4f, .4f, 1.0f));
+        }
+        stateset->setAttributeAndModes(alphaFunc.get());
+        stateset->setAttribute(program.get());
+        stateset->addUniform(baseTextureSampler.get());
+        stateset->setMode(GL_VERTEX_PROGRAM_TWO_SIDE, StateAttribute::ON);
+        stateset->setAttribute(material.get());
+        // XXX This should really come from a material definition
+        // instead of being hard-coded.
+        treeTextureMap.insert(StateSetMap::value_type(firstTree.texture,
+                                                      stateset));
+    } else {
+        stateset = iter->second.get();
+    }
     // Now, create a quadtree for the forest.
-    osg::ref_ptr<osg::Group> _root;
-    ShaderGeometry* leaves[SG_TREE_QUAD_TREE_SIZE][SG_TREE_QUAD_TREE_SIZE];
-
-    // Determine the extents of the tree, and a list of the required textures for later.    
-    osg::BoundingBox extents;
-    for (unsigned int i = 0; i < forest.getNumTrees(); i++)
-    {    
-        const osg::Vec3f center = forest.getTree(i).position.osg() * transform;
-        extents.expandBy(center);
+    {
+        ShaderGeometryQuadtree quadtree(GetTreeCoord(Matrix::inverse(transform)),
+                                        AddTreesLeafObject(),
+                                        SG_TREE_QUAD_TREE_DEPTH,
+                                        MakeTreesLeaf(firstTree.range,
+                                                      shared_geometry));
+        quadtree.buildQuadTree(forest._trees.begin(), forest._trees.end());
+        group = quadtree.getRoot();
     }
-    
-    const osg::Vec2 quadMin(extents.xMin(), extents.yMin());
-    const osg::Vec2 quadMax(extents.xMax(), extents.yMax());
-
-    for (int i = 0; i < SG_TREE_QUAD_TREE_SIZE; ++i) {
-      osg::LOD* interior = new osg::LOD;      
-      //osg::Group* interior = new osg::Group;
-      group->addChild(interior);
-      for (int j = 0; j < SG_TREE_QUAD_TREE_SIZE; ++j) {
-        osg::Geode* geode = new osg::Geode;
-        leaves[i][j] = new ShaderGeometry();
-        leaves[i][j]->setGeometry(shared_geometry);
-        geode->setStateSet(stateset);
-        geode->addDrawable(leaves[i][j]);
-        interior->addChild(geode, 0, firstTree.range);        
-      }
-    }
-    
-    // Now we've got our quadtree, add the trees based on location.
-    
-    for (unsigned int i = 0; i < forest.getNumTrees(); i++)
-    {    
-      TreeBin::Tree t = forest.getTree(i);
-      osg::Vec3 center = t.position.osg() * transform;
-
-      int x = (int)(SG_TREE_QUAD_TREE_SIZE * (center.x() - quadMin.x()) / (quadMax.x() - quadMin.x()));
-      x = osg::clampTo(x, 0, (SG_TREE_QUAD_TREE_SIZE - 1));
-      int y = (int)(SG_TREE_QUAD_TREE_SIZE * (center.y() - quadMin.y()) / (quadMax.y() - quadMin.y()));
-      y = osg::clampTo(y, 0, (SG_TREE_QUAD_TREE_SIZE -1));
-    
-      leaves[y][x]->addTree(t.position.osg(), t.height);
-    }
-    
-    return group;
+    group->setStateSet(stateset);
+    return group.release();    
 }
 
 }
diff --git a/simgear/scene/tgdb/TreeBin.hxx b/simgear/scene/tgdb/TreeBin.hxx
index 48a327fe..da94313e 100644
--- a/simgear/scene/tgdb/TreeBin.hxx
+++ b/simgear/scene/tgdb/TreeBin.hxx
@@ -55,12 +55,10 @@ public:
     { return _trees.size(); }
     const Tree& getTree(unsigned i) const
     { return _trees[i]; }
-
-private:
     TreeList _trees;
 };
 
 osg::Geometry* createOrthQuads(float w, float h, const osg::Matrix& rotate);
-osg::Group* createForest(const TreeBin& forest, const osg::Matrix& transform);
+osg::Group* createForest(TreeBin& forest, const osg::Matrix& transform);
 }
 #endif
diff --git a/simgear/scene/tgdb/obj.cxx b/simgear/scene/tgdb/obj.cxx
index 5fd630d0..d319ec4c 100644
--- a/simgear/scene/tgdb/obj.cxx
+++ b/simgear/scene/tgdb/obj.cxx
@@ -474,7 +474,7 @@ struct SGTileGeometryBin {
       i->second.addRandomSurfacePoints(coverage, 0, randomPoints);
       std::vector<SGVec3f>::iterator j;
       for (j = randomPoints.begin(); j != randomPoints.end(); ++j) {      
-        int k = mt_rand(&seed) * textures.size();        
+        int k = (int)(mt_rand(&seed) * textures.size());        
         if (k == textures.size()) k--;         
         randomForest.insert(*j, textures[k], mat->get_tree_height(), mat->get_tree_width(), mat->get_tree_range());
       }
@@ -515,7 +515,7 @@ struct SGTileGeometryBin {
               i->second.addRandomPoints(object->get_coverage_m2(), randomPoints);
               std::vector<SGVec3f>::iterator l;
               for (l = randomPoints.begin(); l != randomPoints.end(); ++l) {
-                randomModels.insert(*l, object, object->get_randomized_range_m(&seed));
+                randomModels.insert(*l, object, (int)object->get_randomized_range_m(&seed));
               }
             }
           }
@@ -534,6 +534,30 @@ struct SGTileGeometryBin {
   }
 };
 
+typedef std::pair<osg::Node*, int> ModelLOD;
+struct MakeQuadLeaf {
+    osg::LOD* operator() () const { return new osg::LOD; }
+};
+struct AddModelLOD {
+    void operator() (osg::LOD* leaf, ModelLOD& mlod) const
+    {
+        leaf->addChild(mlod.first, 0, mlod.second);
+    }
+};
+struct GetModelLODCoord {
+    GetModelLODCoord(const osg::Matrix& transform) : _transform(transform) {}
+    GetModelLODCoord(const GetModelLODCoord& rhs) : _transform(rhs._transform)
+    {}
+    osg::Vec3 operator() (const ModelLOD& mlod) const
+    {
+        return mlod.first->getBound().center() * _transform;
+    }
+    osg::Matrix _transform;
+};
+
+typedef QuadTreeBuilder<osg::LOD*, ModelLOD, MakeQuadLeaf, AddModelLOD,
+                        GetModelLODCoord>  RandomObjectsQuadtree;
+
 osg::Node*
 SGLoadBTG(const std::string& path, SGMaterialLib *matlib, bool calc_lights, bool use_random_objects)
 {
@@ -549,7 +573,6 @@ SGLoadBTG(const std::string& path, SGMaterialLib *matlib, bool calc_lights, bool
   SGGeod geodPos = SGGeod::fromCart(center);
   SGQuatd hlOr = SGQuatd::fromLonLat(geodPos);
   SGVec3f up = toVec3f(hlOr.backTransform(SGVec3d(0, 0, -1)));
-  osg::Matrix world2Tile(-hlOr.osg());
   GroundLightManager* lightManager = GroundLightManager::instance();
 
   osg::ref_ptr<osg::Group> lightGroup = new SGOffsetTransform(0.94);
@@ -569,6 +592,14 @@ SGLoadBTG(const std::string& path, SGMaterialLib *matlib, bool calc_lights, bool
                                   0, -1,  0, 0,
                                   0,  0, -1, 0,
                                   0,  0,  0, 1);     
+    // Determine an rotation matrix for the models to place them 
+    // perpendicular to the earth's surface. We use the same matrix, 
+    // based on the centre of the tile, as the small angular differences  
+    // between different points on the tile aren't worth worrying about 
+    // for random objects. We also need to flip the orientation 180 degrees
+    osg::Matrix mAtt = flip * osg::Matrix::rotate(hlOr.osg());
+    // The inverse goes from world coordinates to Z up tile coordinates.
+    osg::Matrix world2Tile(osg::Matrix(hlOr.osg().conj()) * flip);
   
     tileGeometryBin.computeRandomObjects(matlib);
     
@@ -576,16 +607,8 @@ SGLoadBTG(const std::string& path, SGMaterialLib *matlib, bool calc_lights, bool
       // Generate a repeatable random seed
       mt seed;
       mt_init(&seed, unsigned(123));
-    
-      // Determine an rotation matrix for the models to place them 
-      // perpendicular to the earth's surface. We use the same matrix, 
-      // based on the centre of the tile, as the small angular differences  
-      // between different points on the tile aren't worth worrying about 
-      // for random objects. We also need to flip the orientation 180 degrees
-      osg::Matrix mAtt = flip * osg::Matrix::rotate(hlOr.osg());
-      
-      LodMap models;
-      
+
+      std::vector<ModelLOD> models;
       for (unsigned int i = 0; i < tileGeometryBin.randomModels.getNumModels(); i++) {
         SGMatModelBin::MatModel obj = tileGeometryBin.randomModels.getMatModel(i);
         osg::Node* node = sgGetRandomModel(obj.model);
@@ -609,10 +632,12 @@ SGLoadBTG(const std::string& path, SGMaterialLib *matlib, bool calc_lights, bool
         }
 
         position->addChild(node);        
-        models.insert(std::pair<osg::ref_ptr<osg::Node>,int>(position, obj.lod));
+        models.push_back(ModelLOD(position, obj.lod));
       }
-      
-      randomObjects = QuadTreeBuilder::makeQuadTree(models, world2Tile);
+      RandomObjectsQuadtree quadtree((GetModelLODCoord(world2Tile)),
+                                     (AddModelLOD()));
+      quadtree.buildQuadTree(models.begin(), models.end());
+      randomObjects = quadtree.getRoot();
       randomObjects->setName("random objects");
     }
       
@@ -620,8 +645,7 @@ SGLoadBTG(const std::string& path, SGMaterialLib *matlib, bool calc_lights, bool
     tileGeometryBin.computeRandomForest(matlib);
     
     if (tileGeometryBin.randomForest.getNumTrees() > 0) {
-      osg::Matrix forAtt = flip * world2Tile;
-      randomForest = createForest(tileGeometryBin.randomForest, forAtt);
+      randomForest = createForest(tileGeometryBin.randomForest, mAtt);
       randomForest->setName("random trees");      
     } 
   }
diff --git a/simgear/scene/util/QuadTreeBuilder.cxx b/simgear/scene/util/QuadTreeBuilder.cxx
index b4dcf072..199efb7d 100644
--- a/simgear/scene/util/QuadTreeBuilder.cxx
+++ b/simgear/scene/util/QuadTreeBuilder.cxx
@@ -25,30 +25,48 @@ using namespace osg;
 
 namespace simgear
 {
-QuadTreeBuilder::QuadTreeBuilder(const Vec2& min, const Vec2& max) :
-    _root(new osg::Group), _min(min), _max(max)
+#if 0
+QuadTreeBuilder::QuadTreeBuilder(const Vec2& min, const Vec2& max, int depth) :
+    _root(new osg::Group), _min(min), _max(max), _depth(depth),
+    _dimension(1 << depth), _leafStorage(_dimension * _dimension),
+    _leaves(_leafStorage, _dimension)
 {
-    for (int i = 0; i < QUAD_TREE_LEAVES; ++i) {
-        Group* interior = new osg::Group;
-        _root->addChild(interior);
-        for (int j = 0; j < QUAD_TREE_LEAVES; ++j) {
-            LOD* lod  = new osg::LOD;
-            interior->addChild(lod);
-            _leaves[i][j] = lod;
+    for (LeafVector::iterator iter = _leafStorage.begin();
+         iter != _leafStorage.end();
+         ++iter)
+        *iter = new LOD;
+    vector<Group*> parentNodes(1);
+    parentNodes[0] = _root.get();
+    unsigned leafDim = 2;
+    for (int i = 0; i < depth - 1;  ++i, leafDim *= 2) {
+        VectorArrayAdapter<vector<Group*> > parents(parentNodes, leafDim / 2);
+        vector<Group*> interiorNodes(leafDim * leafDim);
+        VectorArrayAdapter<vector<Group*> > interiors(interiorNodes, leafDim);
+        for (unsigned j = 0; j < leafDim; ++j) {
+            for (unsigned k = 0; k < leafDim; ++k) {
+                interiors(j, k) = new Group;
+                parents(j / 2, k / 2)->addChild(interiors(j, k));
+            }
         }
+        parentNodes.swap(interiorNodes);
     }
+    VectorArrayAdapter<vector<Group*> > leafParents(parentNodes,
+                                                    _dimension / 2);
+    for (int j = 0; j < _dimension; ++j)
+        for (int k =0; k < _dimension; ++k)
+            leafParents(j / 2, k / 2)->addChild(_leaves(j, k));
 }
 
 void QuadTreeBuilder::addNode(Node* node, int lod, const Matrix& transform)
 {
     Vec3 center = node->getBound().center() * transform;
 
-    int x = (int)(QUAD_TREE_LEAVES * (center.x() - _min.x()) / (_max.x() - _min.x()));
-    x = clampTo(x, 0, (QUAD_TREE_LEAVES - 1));
-    int y = (int)(QUAD_TREE_LEAVES * (center.y() - _min.y()) / (_max.y() - _min.y()));
-    y = clampTo(y, 0, (QUAD_TREE_LEAVES -1));
+    int x = (int)(_dimension * (center.x() - _min.x()) / (_max.x() - _min.x()));
+    x = clampTo(x, 0, (_dimension - 1));
+    int y = (int)(_dimension * (center.y() - _min.y()) / (_max.y() - _min.y()));
+    y = clampTo(y, 0, (_dimension -1));
     
-    _leaves[y][x]->addChild(node, 0, lod);
+    _leaves(y, x)->addChild(node, 0, lod);
 }
 
 osg::Group* QuadTreeBuilder::makeQuadTree(LodMap& models,
@@ -74,5 +92,5 @@ osg::Group* QuadTreeBuilder::makeQuadTree(LodMap& models,
     }
     return result.release();
 }
-
+#endif
 }
diff --git a/simgear/scene/util/QuadTreeBuilder.hxx b/simgear/scene/util/QuadTreeBuilder.hxx
index ff16fb46..2c5c4bf5 100644
--- a/simgear/scene/util/QuadTreeBuilder.hxx
+++ b/simgear/scene/util/QuadTreeBuilder.hxx
@@ -16,36 +16,115 @@
 #ifndef SIMGEAR_QUADTREEBUILDER_HXX
 #define SIMGEAR_QUADTREEBUILDER_HXX 1
 
+#include <algorithm>
+#include <functional>
 #include <vector>
+#include <osg/BoundingBox>
 #include <osg/ref_ptr>
 #include <osg/Node>
 #include <osg/Group>
 #include <osg/Matrix>
 #include <osg/Vec2>
 #include <osg/LOD>
-
-#define QUAD_TREE_LEAVES 4
+#include "VectorArrayAdapter.hxx"
 
 namespace simgear
 {
 typedef std::map<osg::ref_ptr<osg::Node>,int> LodMap;
 
 // Create a quad tree based on x, y extents
+template <typename LeafType, typename ObjectType, typename MakeLeaf,
+          typename AddLeafObject, typename GetObjectLocalCoords>
 class QuadTreeBuilder {
 public:
-    QuadTreeBuilder(const osg::Vec2& min, const osg::Vec2& max);
+    QuadTreeBuilder(const GetObjectLocalCoords& getLocalCoords,
+                    const AddLeafObject& addLeafObject, int depth = 2,
+                    const MakeLeaf& makeLeaf = MakeLeaf()) :
+        _root(new osg::Group), _depth(depth),
+        _dimension(1 << depth), _leafStorage(_dimension * _dimension),
+        _leaves(_leafStorage, _dimension), _getLocalCoords(getLocalCoords),
+        _addLeafObject(addLeafObject), _makeLeaf(makeLeaf)
+    {
+        using namespace std;
+        using namespace osg;
+        generate(_leafStorage.begin(), _leafStorage.end(), _makeLeaf);
+        vector<Group*> parentNodes(1);
+        parentNodes[0] = _root.get();
+        unsigned leafDim = 2;
+        for (int i = 0; i < depth - 1;  ++i, leafDim *= 2) {
+            VectorArrayAdapter<vector<Group*> > parents(parentNodes, leafDim / 2);
+            vector<Group*> interiorNodes(leafDim * leafDim);
+            VectorArrayAdapter<vector<Group*> > interiors(interiorNodes, leafDim);
+            for (unsigned j = 0; j < leafDim; ++j) {
+                for (unsigned k = 0; k < leafDim; ++k) {
+                    interiors(j, k) = new Group;
+                    parents(j / 2, k / 2)->addChild(interiors(j, k));
+                }
+            }
+            parentNodes.swap(interiorNodes);
+        }
+        VectorArrayAdapter<vector<Group*> > leafParents(parentNodes,
+                                                        _dimension / 2);
+        for (int j = 0; j < _dimension; ++j)
+            for (int k =0; k < _dimension; ++k)
+                leafParents(j / 2, k / 2)->addChild(_leaves(j, k));
+    }
+    osg::Vec2 getMin() { return _min; }
+    void setMin(const osg::Vec2& min) { _min = min; }
+    osg::Vec2 getMax() { return _max; }
+    void setMax(const osg::Vec2& max) { _max = max; }
     ~QuadTreeBuilder() {}
     osg::Group* getRoot() { return _root.get(); }
-    // Add node to the quadtree using its x, y and LoD
-    void addNode(osg::Node* node, int lod, const osg::Matrix& transform);
+
+    void addNode(ObjectType& obj)
+    {
+        using namespace osg;
+        const Vec3 center(_getLocalCoords(obj));
+        int x = (int)(_dimension * (center.x() - _min.x())
+                      / (_max.x() - _min.x()));
+        x = clampTo(x, 0, (_dimension - 1));
+        int y = (int)(_dimension * (center.y() - _min.y())
+                      / (_max.y() - _min.y()));
+        y = clampTo(y, 0, (_dimension -1));
+        _addLeafObject(_leaves(y, x), obj);
+    }
+    // STL craziness
+    struct AddNode
+    {
+        AddNode(QuadTreeBuilder* qt) : _qt(qt) {}
+        AddNode(const AddNode& rhs) : _qt(rhs._qt) {}
+        void operator() (ObjectType& obj) const { _qt->addNode(obj); }
+        QuadTreeBuilder *_qt;
+    };
     // Make a quadtree of nodes from a map of nodes and LOD values
-    static osg::Group* makeQuadTree(LodMap& nodes,
-                                    const osg::Matrix& transform);
+    template <typename ForwardIterator>
+    void buildQuadTree(const ForwardIterator& begin,
+                       const ForwardIterator& end)
+    {
+        using namespace osg;
+        BoundingBox extents;
+        for (ForwardIterator iter = begin; iter != end; ++iter) {
+            const Vec3 center = _getLocalCoords(*iter);
+            extents.expandBy(center);
+        }
+        _min = Vec2(extents.xMin(), extents.yMin());
+        _max = Vec2(extents.xMax(), extents.yMax());
+        std::for_each(begin, end, AddNode(this));
+    }
+
 protected:
+    typedef std::vector<LeafType> LeafVector;
     osg::ref_ptr<osg::Group> _root;
-    osg::LOD* _leaves[QUAD_TREE_LEAVES][QUAD_TREE_LEAVES];
     osg::Vec2 _min;
     osg::Vec2 _max;
+    int _depth;
+    int _dimension;
+    LeafVector _leafStorage;
+    VectorArrayAdapter<LeafVector> _leaves;
+    const GetObjectLocalCoords _getLocalCoords;
+    const AddLeafObject _addLeafObject;
+    const MakeLeaf _makeLeaf;
 };
+
 }
 #endif
diff --git a/simgear/scene/util/RenderConstants.hxx b/simgear/scene/util/RenderConstants.hxx
index 66ce522e..184b9abc 100644
--- a/simgear/scene/util/RenderConstants.hxx
+++ b/simgear/scene/util/RenderConstants.hxx
@@ -50,6 +50,10 @@ enum NodeMask {
 //
 // Normal opaque objects are assigned bin 0.
 //
+// Random objects like trees may have transparency, but there are too
+// many to depth sort. By drawing them after the terrain we can at
+// least keep the sky under the ground from poking through.
+//
 // Point lights blend with the terrain to simulate attenuation but
 // should completely obscure any transparent geometry behind
 // them. Also, they should be visible through semi-transparent cloud
@@ -62,6 +66,7 @@ enum NodeMask {
 // OSG and its file loaders throw all transparent objects into bin 10.
 
 enum RenderBin {
+    RANDOM_OBJECTS_BIN = 2,
     POINT_LIGHTS_BIN = 8,
     CLOUDS_BIN = 9,
     TRANSPARENT_BIN = 10        // assigned by OSG
diff --git a/simgear/scene/util/StateAttributeFactory.cxx b/simgear/scene/util/StateAttributeFactory.cxx
index 17a1e305..30407b05 100644
--- a/simgear/scene/util/StateAttributeFactory.cxx
+++ b/simgear/scene/util/StateAttributeFactory.cxx
@@ -57,6 +57,11 @@ StateAttributeFactory::StateAttributeFactory()
     _whiteTexture->setWrap(osg::Texture::WRAP_S, osg::Texture::REPEAT);
     _whiteTexture->setWrap(osg::Texture::WRAP_T, osg::Texture::REPEAT);
     _whiteTexture->setDataVariance(osg::Object::STATIC);
+    _white = new Vec4Array(1);
+    (*_white)[0].set(1.0f, 1.0f, 1.0f, 1.0f);
+    _white->setDataVariance(Object::STATIC);
+    _cullFaceBack = new CullFace(CullFace::BACK);
+    _cullFaceBack->setDataVariance(Object::STATIC);
 }
 
 osg::ref_ptr<StateAttributeFactory> StateAttributeFactory::_theInstance;
diff --git a/simgear/scene/util/StateAttributeFactory.hxx b/simgear/scene/util/StateAttributeFactory.hxx
index 080a386f..c638108c 100644
--- a/simgear/scene/util/StateAttributeFactory.hxx
+++ b/simgear/scene/util/StateAttributeFactory.hxx
@@ -25,7 +25,9 @@
 #include <OpenThreads/Mutex>
 #include <osg/ref_ptr>
 #include <osg/AlphaFunc>
+#include <osg/Array>
 #include <osg/BlendFunc>
+#include <osg/CullFace>
 #include <osg/ShadeModel>
 #include <osg/Texture2D>
 #include <osg/TexEnv>
@@ -45,6 +47,11 @@ public:
     osg::ShadeModel* getFlatShadeModel() { return _flat.get(); }
     // White, repeating texture
     osg::Texture2D* getWhiteTexture() { return _whiteTexture.get(); }
+    // White color
+    osg::Vec4Array* getWhiteColor() {return _white.get(); }
+    // cull back facing polygons
+    osg::CullFace* getCullFaceBack() { return _cullFaceBack.get(); }
+    
     static StateAttributeFactory* instance();
 protected:
     StateAttributeFactory();
@@ -54,6 +61,8 @@ protected:
     osg::ref_ptr<osg::BlendFunc> _standardBlendFunc;
     osg::ref_ptr<osg::TexEnv> _standardTexEnv;
     osg::ref_ptr<osg::Texture2D> _whiteTexture;
+    osg::ref_ptr<osg::Vec4Array> _white;
+    osg::ref_ptr<osg::CullFace> _cullFaceBack;
     static osg::ref_ptr<StateAttributeFactory> _theInstance;
     static OpenThreads::Mutex _instanceMutex;
 };
diff --git a/simgear/scene/util/VectorArrayAdapter.hxx b/simgear/scene/util/VectorArrayAdapter.hxx
index 44c515b9..700834b2 100644
--- a/simgear/scene/util/VectorArrayAdapter.hxx
+++ b/simgear/scene/util/VectorArrayAdapter.hxx
@@ -22,6 +22,7 @@
 #ifndef VECTORARRAYADAPTERHXX
 #define VECTORARRAYADAPTERHXX 1
 
+// #define SG_CHECK_VECTOR_ACCESS 1
 namespace simgear
 {
 template <typename Vector>
@@ -44,7 +45,16 @@ public:
         _rowOffset(rowOffset)
     {
     }
-    
+#ifdef SG_CHECK_VECTOR_ACCESS
+    typename Vector::value_type& operator() (int i, int j)
+    {
+        return _v.at(_baseOffset + i * _rowStride + _rowOffset + j);
+    }
+    const typename Vector::value_type& operator() (int i, int j) const
+    {
+        return _v.at(_baseOffset + i * _rowStride + _rowOffset + j);
+    }
+#else
     typename Vector::value_type& operator() (int i, int j)
     {
         return _v[_baseOffset + i * _rowStride + _rowOffset + j];
@@ -53,6 +63,7 @@ public:
     {
         return _v[_baseOffset + i * _rowStride + _rowOffset + j];
     }
+#endif
 private:
     Vector& _v;
     const int _rowStride;