diff --git a/src/assembly.rs b/src/assembly.rs
index 44484e5..7d8d659 100644
--- a/src/assembly.rs
+++ b/src/assembly.rs
@@ -9,6 +9,7 @@ use light_accel::{LightAccel, LightTree};
 use light::LightSource;
 use math::{Matrix4x4, Vector};
 use surface::{Surface, SurfaceIntersection};
+use transform_stack::TransformStack;
 
 
 #[derive(Debug)]
@@ -32,30 +33,56 @@ pub struct Assembly {
 }
 
 impl Assembly {
-    // Returns (light_color, shadow_vector, selection_pdf)
+    // Returns (light_color, shadow_vector, pdf, selection_pdf)
     pub fn sample_lights(&self,
+                         xform_stack: &mut TransformStack,
                          n: f32,
                          uvw: (f32, f32, f32),
                          wavelength: f32,
                          time: f32,
                          intr: &SurfaceIntersection)
-                         -> Option<(SpectralSample, Vector, f32)> {
+                         -> Option<(SpectralSample, Vector, f32, f32)> {
         if let &SurfaceIntersection::Hit { pos, incoming, nor, closure, .. } = intr {
-            if let Some((light_i, sel_pdf, _)) = self.light_accel
-                .select(incoming, pos, nor, closure.as_surface_closure(), time, n) {
+            let sel_xform = if xform_stack.top().len() > 0 {
+                lerp_slice(xform_stack.top(), time)
+            } else {
+                Matrix4x4::new()
+            };
+            if let Some((light_i, sel_pdf, whittled_n)) = self.light_accel
+                .select(incoming * sel_xform,
+                        pos * sel_xform,
+                        nor * sel_xform,
+                        closure.as_surface_closure(),
+                        time,
+                        n) {
                 let inst = self.light_instances[light_i];
                 match inst.instance_type {
+
                     InstanceType::Object => {
                         match &self.objects[inst.data_index] {
                             &Object::Light(ref light) => {
+                                // Get the world-to-object space transform of the light
                                 let xform = if let Some((a, b)) = inst.transform_indices {
-                                    lerp_slice(&self.xforms[a..b], time)
+                                    let pxforms = xform_stack.top();
+                                    let xform = lerp_slice(&self.xforms[a..b], time);
+                                    if pxforms.len() > 0 {
+                                        lerp_slice(pxforms, time) * xform
+                                    } else {
+                                        xform
+                                    }
                                 } else {
-                                    Matrix4x4::new()
+                                    let pxforms = xform_stack.top();
+                                    if pxforms.len() > 0 {
+                                        lerp_slice(pxforms, time)
+                                    } else {
+                                        Matrix4x4::new()
+                                    }
                                 };
+
+                                // Sample the light
                                 let (color, shadow_vec, pdf) =
                                     light.sample(&xform, pos, uvw.0, uvw.1, wavelength, time);
-                                return Some((color, shadow_vec, pdf * sel_pdf));
+                                return Some((color, shadow_vec, pdf, sel_pdf));
                             }
 
                             _ => unimplemented!(),
@@ -63,9 +90,23 @@ impl Assembly {
                     }
 
                     InstanceType::Assembly => {
-                        // TODO: recursive light selection inside assemblies
-                        unimplemented!()
+                        // Push the world-to-object space transforms of the assembly onto
+                        // the transform stack.
+                        if let Some((a, b)) = inst.transform_indices {
+                            xform_stack.push(&self.xforms[a..b]);
+                        }
 
+                        // Sample sub-assembly lights
+                        let sample = self.assemblies[inst.data_index]
+                            .sample_lights(xform_stack, whittled_n, uvw, wavelength, time, intr);
+
+                        // Pop the assembly's transforms off the transform stack.
+                        if let Some(_) = inst.transform_indices {
+                            xform_stack.pop();
+                        }
+
+                        // Return sample
+                        return sample.map(|(ss, v, pdf, spdf)| (ss, v, pdf, spdf * sel_pdf));
                     }
                 }
             } else {
@@ -190,8 +231,7 @@ impl AssemblyBuilder {
         self.objects.shrink_to_fit();
         self.assemblies.shrink_to_fit();
 
-        // Calculate instance bounds, used for building object accel and
-        // (TODO) light accel.
+        // Calculate instance bounds, used for building object accel and light accel.
         let (bis, bbs) = self.instance_bounds();
 
         // Build object accel
@@ -199,8 +239,8 @@ impl AssemblyBuilder {
                                              1,
                                              |inst| &bbs[bis[inst.id]..bis[inst.id + 1]]);
 
-        // Get list of instances that are for light sources.
-        // TODO: include assemblies that themselves contain light sources.
+        // Get list of instances that are for light sources or assemblies that contain light
+        // sources.
         let mut light_instances: Vec<_> = self.instances
             .iter()
             .filter(|inst| {
@@ -213,7 +253,9 @@ impl AssemblyBuilder {
                         }
                     }
 
-                    _ => false,
+                    InstanceType::Assembly => {
+                        self.assemblies[inst.data_index].light_accel.approximate_energy() > 0.0
+                    }
                 }
             })
             .map(|&a| a)
@@ -231,8 +273,9 @@ impl AssemblyBuilder {
                     }
                 }
 
-                // TODO: handle assemblies.
-                _ => 0.0,
+                InstanceType::Assembly => {
+                    self.assemblies[inst.data_index].light_accel.approximate_energy()
+                }
             };
             (bounds, energy)
         });
diff --git a/src/bvh.rs b/src/bvh.rs
index 2427380..dad822a 100644
--- a/src/bvh.rs
+++ b/src/bvh.rs
@@ -1,15 +1,12 @@
 #![allow(dead_code)]
 
-use std;
-use std::cmp::Ordering;
-
-use algorithm::{partition, quick_select, merge_slices_append};
+use algorithm::{partition, merge_slices_append};
 use bbox::BBox;
 use boundable::Boundable;
 use lerp::lerp_slice;
 use math::log2_64;
 use ray::AccelRay;
-use sah::sah_split;
+use objects_split::{sah_split, median_split};
 
 
 const BVH_MAX_DEPTH: usize = 64;
@@ -125,15 +122,6 @@ impl BVH {
                 split_axis: 0,
             });
 
-            // Get combined object bounds
-            let bounds = {
-                let mut bb = BBox::new();
-                for obj in &objects[..] {
-                    bb |= lerp_slice(bounder(obj), 0.5);
-                }
-                bb
-            };
-
             // Partition objects.
             // If we're too near the max depth, we do balanced building to
             // avoid exceeding max depth.
@@ -144,45 +132,7 @@ impl BVH {
                 sah_split(objects, &bounder)
             } else {
                 // Balanced splitting, when we don't have room to play
-                let split_axis = {
-                    let mut axis = 0;
-                    let mut largest = std::f32::NEG_INFINITY;
-                    for i in 0..3 {
-                        let extent = bounds.max.get_n(i) - bounds.min.get_n(i);
-                        if extent > largest {
-                            largest = extent;
-                            axis = i;
-                        }
-                    }
-                    axis
-                };
-
-                let place = {
-                    let place = objects.len() / 2;
-                    if place > 0 {
-                        place
-                    } else {
-                        1
-                    }
-                };
-                quick_select(objects, place, |a, b| {
-                    let tb_a = lerp_slice(bounder(a), 0.5);
-                    let tb_b = lerp_slice(bounder(b), 0.5);
-                    let centroid_a = (tb_a.min.get_n(split_axis) + tb_a.max.get_n(split_axis)) *
-                                     0.5;
-                    let centroid_b = (tb_b.min.get_n(split_axis) + tb_b.max.get_n(split_axis)) *
-                                     0.5;
-
-                    if centroid_a < centroid_b {
-                        Ordering::Less
-                    } else if centroid_a == centroid_b {
-                        Ordering::Equal
-                    } else {
-                        Ordering::Greater
-                    }
-                });
-
-                (place, split_axis)
+                median_split(objects, &bounder)
             };
 
             // Create child nodes
diff --git a/src/light/sphere_light.rs b/src/light/sphere_light.rs
index 25ba4a3..4f43308 100644
--- a/src/light/sphere_light.rs
+++ b/src/light/sphere_light.rs
@@ -44,8 +44,10 @@ impl LightSource for SphereLight {
               wavelength: f32,
               time: f32)
               -> (SpectralSample, Vector, f32) {
-        // TODO: use transform space correctly
-        let pos = Point::new(0.0, 0.0, 0.0) * space.inverse();
+        // TODO: track fp error due to transforms
+        let arr = arr * *space;
+        let pos = Point::new(0.0, 0.0, 0.0);
+
         // Calculate time interpolated values
         let radius: f64 = lerp_slice(&self.radii, time) as f64;
         let col = lerp_slice(&self.colors, time);
@@ -93,13 +95,14 @@ impl LightSource for SphereLight {
                                      (d - (cos_a * radius)) as f32);
 
             // Calculate the final values and return everything.
-            let shadow_vec = (x * sample.x()) + (y * sample.y()) + (z * sample.z());
+            let shadow_vec = ((x * sample.x()) + (y * sample.y()) + (z * sample.z())) *
+                             space.inverse();
             let pdf = uniform_sample_cone_pdf(cos_theta_max);
             let spectral_sample = (col * surface_area_inv as f32).to_spectral_sample(wavelength);
             return (spectral_sample, shadow_vec, pdf as f32);
         } else {
             // If we're inside the sphere, there's light from every direction.
-            let shadow_vec = uniform_sample_sphere(u, v);
+            let shadow_vec = uniform_sample_sphere(u, v) * space.inverse();
             let pdf = 1.0 / (4.0 * PI_64);
             let spectral_sample = (col * surface_area_inv as f32).to_spectral_sample(wavelength);
             return (spectral_sample, shadow_vec, pdf as f32);
@@ -118,8 +121,8 @@ impl LightSource for SphereLight {
         // We're not using these, silence warnings
         let _ = (sample_dir, sample_u, sample_v, wavelength);
 
-        // TODO: use transform space correctly
-        let pos = Point::new(0.0, 0.0, 0.0) * space.inverse();
+        let arr = arr * *space;
+        let pos = Point::new(0.0, 0.0, 0.0);
         let radius: f64 = lerp_slice(&self.radii, time) as f64;
 
         let d2: f64 = (pos - arr).length2() as f64;  // Distance from center of sphere squared
diff --git a/src/light_accel/light_tree.rs b/src/light_accel/light_tree.rs
index 2f9896e..97dd2c1 100644
--- a/src/light_accel/light_tree.rs
+++ b/src/light_accel/light_tree.rs
@@ -2,7 +2,7 @@ use algorithm::merge_slices_append;
 use bbox::BBox;
 use lerp::lerp_slice;
 use math::{Vector, Point, Normal};
-use sah::sah_split;
+use objects_split::sah_split;
 use shading::surface_closure::SurfaceClosure;
 
 use super::LightAccel;
@@ -187,4 +187,12 @@ impl LightAccel for LightTree {
         // Found our light!
         Some((self.nodes[node_index].child_index, tot_prob, n))
     }
+
+    fn approximate_energy(&self) -> f32 {
+        if self.nodes.len() > 0 {
+            self.nodes[0].energy
+        } else {
+            0.0
+        }
+    }
 }
diff --git a/src/light_accel/mod.rs b/src/light_accel/mod.rs
index 1dbeba6..7d42f6f 100644
--- a/src/light_accel/mod.rs
+++ b/src/light_accel/mod.rs
@@ -17,11 +17,14 @@ pub trait LightAccel {
               time: f32,
               n: f32)
               -> Option<(usize, f32, f32)>;
+
+    fn approximate_energy(&self) -> f32;
 }
 
 #[derive(Debug, Clone)]
 pub struct LightArray {
     indices: Vec<usize>,
+    aprx_energy: f32,
 }
 
 impl LightArray {
@@ -30,15 +33,20 @@ impl LightArray {
         where F: 'a + Fn(&T) -> Option<(&'a [BBox], f32)>
     {
         let mut indices = Vec::new();
+        let mut aprx_energy = 0.0;
         for (i, thing) in things.iter().enumerate() {
             if let Some((_, power)) = q(thing) {
                 if power > 0.0 {
                     indices.push(i);
+                    aprx_energy += power;
                 }
             }
         }
 
-        LightArray { indices: indices }
+        LightArray {
+            indices: indices,
+            aprx_energy: aprx_energy,
+        }
     }
 }
 
@@ -71,4 +79,8 @@ impl LightAccel for LightArray {
 
         Some((i, pdf, whittled_n))
     }
+
+    fn approximate_energy(&self) -> f32 {
+        self.aprx_energy
+    }
 }
diff --git a/src/main.rs b/src/main.rs
index 7f0e2f4..e126247 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -28,10 +28,10 @@ mod lerp;
 mod light_accel;
 mod light;
 mod math;
+mod objects_split;
 mod parse;
 mod ray;
 mod renderer;
-mod sah;
 mod sampling;
 mod scene;
 mod shading;
diff --git a/src/sah.rs b/src/objects_split.rs
similarity index 50%
rename from src/sah.rs
rename to src/objects_split.rs
index 51168cb..b1a4985 100644
--- a/src/sah.rs
+++ b/src/objects_split.rs
@@ -1,6 +1,7 @@
 use std;
+use std::cmp::Ordering;
 
-use algorithm::partition;
+use algorithm::{partition, quick_select};
 use bbox::BBox;
 use lerp::lerp_slice;
 
@@ -95,3 +96,105 @@ pub fn sah_split<'a, T, F>(objects: &mut [T], bounder: &F) -> (usize, usize)
 
     (split_i, split_axis)
 }
+
+/// Takes a slice of boundable objects and partitions them based on the bounds mean heuristic.
+///
+/// Returns the index of the partition boundary and the axis that it split on
+/// (0 = x, 1 = y, 2 = z).
+pub fn bounds_mean_split<'a, T, F>(objects: &mut [T], bounder: &F) -> (usize, usize)
+    where F: Fn(&T) -> &'a [BBox]
+{
+    // Get combined object bounds
+    let bounds = {
+        let mut bb = BBox::new();
+        for obj in &objects[..] {
+            bb |= lerp_slice(bounder(obj), 0.5);
+        }
+        bb
+    };
+
+    let split_axis = {
+        let mut axis = 0;
+        let mut largest = std::f32::NEG_INFINITY;
+        for i in 0..3 {
+            let extent = bounds.max.get_n(i) - bounds.min.get_n(i);
+            if extent > largest {
+                largest = extent;
+                axis = i;
+            }
+        }
+        axis
+    };
+
+    let div = (bounds.min.get_n(split_axis) + bounds.max.get_n(split_axis)) * 0.5;
+
+    // Partition
+    let mut split_i = partition(&mut objects[..], |obj| {
+        let tb = lerp_slice(bounder(obj), 0.5);
+        let centroid = (tb.min.get_n(split_axis) + tb.max.get_n(split_axis)) * 0.5;
+        centroid < div
+    });
+    if split_i < 1 {
+        split_i = 1;
+    } else if split_i >= objects.len() {
+        split_i = objects.len() - 1;
+    }
+
+    (split_i, split_axis)
+}
+
+
+/// Takes a slice of boundable objects and partitions them based on the median heuristic.
+///
+/// Returns the index of the partition boundary and the axis that it split on
+/// (0 = x, 1 = y, 2 = z).
+pub fn median_split<'a, T, F>(objects: &mut [T], bounder: &F) -> (usize, usize)
+    where F: Fn(&T) -> &'a [BBox]
+{
+    // Get combined object bounds
+    let bounds = {
+        let mut bb = BBox::new();
+        for obj in &objects[..] {
+            bb |= lerp_slice(bounder(obj), 0.5);
+        }
+        bb
+    };
+
+    let split_axis = {
+        let mut axis = 0;
+        let mut largest = std::f32::NEG_INFINITY;
+        for i in 0..3 {
+            let extent = bounds.max.get_n(i) - bounds.min.get_n(i);
+            if extent > largest {
+                largest = extent;
+                axis = i;
+            }
+        }
+        axis
+    };
+
+    let place = {
+        let place = objects.len() / 2;
+        if place > 0 {
+            place
+        } else {
+            1
+        }
+    };
+    quick_select(objects, place, |a, b| {
+        let tb_a = lerp_slice(bounder(a), 0.5);
+        let tb_b = lerp_slice(bounder(b), 0.5);
+        let centroid_a = (tb_a.min.get_n(split_axis) + tb_a.max.get_n(split_axis)) * 0.5;
+        let centroid_b = (tb_b.min.get_n(split_axis) + tb_b.max.get_n(split_axis)) * 0.5;
+
+        if centroid_a < centroid_b {
+            Ordering::Less
+        } else if centroid_a == centroid_b {
+            Ordering::Equal
+        } else {
+            Ordering::Greater
+        }
+    });
+
+    (place, split_axis)
+}
diff --git a/src/renderer.rs b/src/renderer.rs
index 27ac69f..f66f3a9 100644
--- a/src/renderer.rs
+++ b/src/renderer.rs
@@ -18,6 +18,7 @@ use ray::Ray;
 use scene::Scene;
 use surface;
 use tracer::Tracer;
+use transform_stack::TransformStack;
 
 
 #[derive(Debug)]
@@ -67,6 +68,7 @@ impl Renderer {
                     let mut paths = Vec::new();
                     let mut rays = Vec::new();
                     let mut tracer = Tracer::from_assembly(&self.scene.root);
+                    let mut xform_stack = TransformStack::new();
 
                     loop {
                         paths.clear();
@@ -130,7 +132,7 @@ impl Renderer {
                             // Determine next rays to shoot based on result
                             pi =
                                 partition_pair(&mut paths[..pi], &mut rays[..pi], |i, path, ray| {
-                                    path.next(&self.scene, &isects[i], &mut *ray)
+                                    path.next(&mut xform_stack, &self.scene, &isects[i], &mut *ray)
                                 });
                         }
 
@@ -274,7 +276,12 @@ impl LightPath {
         s
     }
 
-    fn next(&mut self, scene: &Scene, isect: &surface::SurfaceIntersection, ray: &mut Ray) -> bool {
+    fn next(&mut self,
+            xform_stack: &mut TransformStack,
+            scene: &Scene,
+            isect: &surface::SurfaceIntersection,
+            ray: &mut Ray)
+            -> bool {
         self.round += 1;
 
         // Result of shading ray, prepare light ray
@@ -291,14 +298,20 @@ impl LightPath {
                 // Prepare light ray
                 let light_n = self.next_lds_samp();
                 let light_uvw = (self.next_lds_samp(), self.next_lds_samp(), self.next_lds_samp());
-                if let Some((light_color, shadow_vec, light_pdf)) = scene.root
-                    .sample_lights(light_n, light_uvw, self.wavelength, self.time, isect) {
+                xform_stack.clear();
+                if let Some((light_color, shadow_vec, light_pdf, light_sel_pdf)) = scene.root
+                    .sample_lights(xform_stack,
+                                   light_n,
+                                   light_uvw,
+                                   self.wavelength,
+                                   self.time,
+                                   isect) {
                     // Calculate and store the light that will be contributed
                     // to the film plane if the light is not in shadow.
                     self.pending_color_addition = {
                         let material = closure.as_surface_closure();
                         let la = material.evaluate(ray.dir, shadow_vec, nor, self.wavelength);
-                        light_color * la * self.light_attenuation / light_pdf
+                        light_color * la * self.light_attenuation / (light_pdf * light_sel_pdf)
                     };
 
                     // Calculate the shadow ray for testing if the light is
diff --git a/src/shading/surface_closure.rs b/src/shading/surface_closure.rs
index d76a5d5..c893a6a 100644
--- a/src/shading/surface_closure.rs
+++ b/src/shading/surface_closure.rs
@@ -323,7 +323,7 @@ impl SurfaceClosure for LambertClosure {
                 }
                 .into_vector();
 
-            let cos_nv = dot(nn, v);
+            let cos_nv = dot(nn, v).max(-1.0).min(1.0);
 
             return sphere_lambert(cos_nv, cos_theta);
         }
diff --git a/todo.txt b/todo.txt
index c93b8f3..64461b3 100644
--- a/todo.txt
+++ b/todo.txt
@@ -1,7 +1,7 @@
 //- Implement support for multiple light sources in a scene.
 //- Implement light tree.
-- Implement proper light source selection for hierarchical instancing.
+//- Implement proper light source selection for hierarchical instancing.
 - Implement shader parsing, and use in rendering.
 - Implement floating point accuracy stuff.
 - Add gui display of rendering in progress.
-- Unit tests for scene parsing.
\ No newline at end of file
+- Unit tests for scene parsing.