sphere works again, initial support for coherent RT with non-tri shapes

SConstruct

@@ -517,7 +517,7 @@ plugins += env.SharedLibrary('plugins/obj', ['src/shapes/obj.cpp'])
 plugins += env.SharedLibrary('plugins/ply', ['src/shapes/ply/ply.cpp', 'src/shapes/ply/ply_parser.cpp'],
 	CPPPATH = env['CPPPATH'] + ['src/shapes/ply'])
 plugins += env.SharedLibrary('plugins/serialized', ['src/shapes/serialized.cpp'])
-#plugins += env.SharedLibrary('plugins/sphere', ['src/shapes/sphere.cpp'])
+plugins += env.SharedLibrary('plugins/sphere', ['src/shapes/sphere.cpp'])
 #plugins += env.SharedLibrary('plugins/cylinder', ['src/shapes/cylinder.cpp'])
 #plugins += env.SharedLibrary('plugins/hair', ['src/shapes/hair.cpp', 'src/shapes/miterseg.cpp'])
 #plugins += env.SharedLibrary('plugins/group', ['src/shapes/group.cpp'])

doc/integrator.tex

@@ -200,7 +200,7 @@ To avoid having to do this every time \code{Li()} is called,
 we can override the \code{preprocess} function:
 \begin{cpp}
 	/// Preprocess function -- called on the initiating machine
-	void preprocess(const Scene *scene, RenderQueue *queue, 
+	bool preprocess(const Scene *scene, RenderQueue *queue, 
 			const RenderJob *job, int sceneResID, int cameraResID, 
 			int samplerResID) {
 		SampleIntegrator::preprocess(scene, queue, job, sceneResID, 
@@ -213,6 +213,8 @@ we can override the \code{preprocess} function:
 		for (int i=0; i<8; ++i)
 			m_maxDist = std::max(m_maxDist, 
 				(cameraPosition - sceneAABB.getCorner(i)).length());
+
+		return true;
 	}
 \end{cpp}
 The bottom of this function should be relatively self-explanatory. The

include/mitsuba/core/ray.h

@@ -185,7 +185,7 @@ struct RayInterval4 {
 		mint = SSEConstants::eps;
 		maxt = SSEConstants::p_inf;
 	}
-	
+
 	inline RayInterval4(const Ray *rays) {
 		for (int i=0; i<4; i++) {
 			mint.f[i] = rays[i].mint;

include/mitsuba/render/gkdtree.h

@@ -40,8 +40,8 @@
 #define MTS_KD_BLOCKSIZE_KD  (512*1024/sizeof(KDNode))
 #define MTS_KD_BLOCKSIZE_IDX (512*1024/sizeof(uint32_t))
 
-/// 64 byte temporary storage for intersection computations 
+/// 32 byte temporary storage for intersection computations 
-#define MTS_KD_INTERSECTION_TEMP 64
+#define MTS_KD_INTERSECTION_TEMP 32
 
 /// Use a simple hashed 8-entry mailbox per thread
 //#define MTS_KD_MAILBOX_ENABLED 1

include/mitsuba/render/kdtree.h

@@ -90,14 +90,14 @@ public:
 	 * use of ray coherence to do this very efficiently. Requires SSE.
 	 */
 	void rayIntersectPacket(const RayPacket4 &packet, 
-		const RayInterval4 &interval, Intersection4 &its) const;
+		const RayInterval4 &interval, Intersection4 &its, void *temp) const;
 
 	/**
 	 * \brief Fallback for incoherent rays
 	 * \sa rayIntesectPacket
 	 */
 	void rayIntersectPacketIncoherent(const RayPacket4 &packet, 
-		const RayInterval4 &interval, Intersection4 &its) const;
+		const RayInterval4 &interval, Intersection4 &its, void *temp) const;
 #endif
 
 	MTS_DECLARE_CLASS()
@@ -135,7 +135,7 @@ protected:
 			const TriMesh *mesh = static_cast<const TriMesh *>(shape);
 			return mesh->getTriangles()[idx].getClippedAABB(mesh->getVertexPositions(), aabb);
 		} else {
-			return shape->getAABB();
+			return shape->getClippedAABB(aabb);
 		}
 	}
 
@@ -146,7 +146,7 @@ protected:
 	/// Temporarily holds some intersection information
 	struct IntersectionCache {
 		size_type shapeIndex;
-		size_type index;
+		size_type primIndex;
 		Float u, v;
 	};
 
@@ -173,7 +173,7 @@ protected:
 					return ENo;
 				t = tempT;
 				cache->shapeIndex = shapeIdx;
-				cache->index = idx;
+				cache->primIndex = primIdx;
 				cache->u = tempU;
 				cache->v = tempV;
 				return EYes;
@@ -181,8 +181,9 @@ protected:
 		} else {
 			const Shape *shape = m_shapes[shapeIndex];
 			if (shape->rayIntersect(ray, mint, maxt, t, 
-						reinterpret_cast<uint8_t*>(temp) + 4)) {
+						reinterpret_cast<uint8_t*>(temp) + 8)) {
 				cache->shapeIndex = shapeIdx;
+				cache->primIndex = KNoTriangleFlag;
 				return EYes;
 			}
 		}
@@ -193,7 +194,7 @@ protected:
 			if (ta.rayIntersect(ray, mint, maxt, tempU, tempV, tempT)) {
 				t = tempT;
 				cache->shapeIndex = ta.shapeIndex;
-				cache->index = ta.index;
+				cache->primIndex = ta.primIndex;
 				cache->u = tempU;
 				cache->v = tempV;
 				return EYes;
@@ -202,8 +203,9 @@ protected:
 			uint32_t shapeIndex = ta.shapeIndex;
 			const Shape *shape = m_shapes[shapeIndex];
 			if (shape->rayIntersect(ray, mint, maxt, t, 
-						reinterpret_cast<uint8_t*>(temp) + 4)) {
+						reinterpret_cast<uint8_t*>(temp) + 8)) {
 				cache->shapeIndex = shapeIndex;
+				cache->primIndex = KNoTriangleFlag;
 				return EYes;
 			}
 		}
@@ -221,17 +223,6 @@ protected:
 	};
 #endif
 
-	/**
-	 * Fallback for incoherent rays
-	 */
-	inline void rayIntersectPacketIncoherent(const Ray *rays, 
-			Intersection *its) const {
-		for (int i=0; i<4; i++) {
-			if (!rayIntersect(rays[i], its[i]))
-				its[i].t = std::numeric_limits<float>::infinity();
-		}
-	}
-
 	/// Virtual destructor
 	virtual ~KDTree();
 private:

include/mitsuba/render/shape.h

@@ -111,9 +111,6 @@ public:
 	/// Return a string representation
 	std::string toString() const;
 public:
-	/// Incident direction in the local frame
-	Vector wi;
-
 	/// Distance traveled along the ray
 	Float t;
 
@@ -129,7 +126,7 @@ public:
 	/// UV surface coordinates
 	Point2 uv;
 
-	/// Position partials wrt. to changes in texture-space
+	/// Position partials wrt. the texture space parameterization
 	Vector dpdu, dpdv;
 
 	/// Texture coordinate mapping partials wrt. changes in screen-space
@@ -138,6 +135,9 @@ public:
 	/// Interpolated vertex color
 	Spectrum color;
 
+	/// Incident direction in the local frame
+	Vector wi;
+
 	/// Affected shape
 	const Shape *shape;
 
@@ -197,13 +197,6 @@ public:
 	virtual void fillIntersectionRecord(const Ray &ray, Float t, 
 			const void *temp, Intersection &its) const;
 
-#if defined(MTS_SSE)
-	/// Perform 4 simultaneous intersection tests using SSE
-//	virtual __m128 rayIntersectPacket(const RayPacket4 &packet, const
-  //      __m128 mint, __m128 maxt, __m128 inactive,
-  //      Intersection4 &its) const;
-#endif
-
 	/**
 	 * \brief Sample a point on the shape
 	 *

include/mitsuba/render/triaccel.h

@@ -44,8 +44,8 @@ struct TriAccel {
 
 	Float c_nu;
 	Float c_nv;
-	uint32_t index;
 	uint32_t shapeIndex;
+	uint32_t primIndex;
 
 	/// Construct from vertex data. Returns '1' if there was a failure
 	inline int load(const Point &A, const Point &B, const Point &C);

src/librender/kdtree.cpp

@@ -86,7 +86,7 @@ void KDTree::build() {
 				const Point &v2 = positions[tri.idx[2]];
 				m_triAccel[idx].load(v0, v1, v2);
 				m_triAccel[idx].shapeIndex = i;
-				m_triAccel[idx].index = j;
+				m_triAccel[idx].primIndex = j;
 				++idx;
 			}
 		} else {
@@ -128,7 +128,7 @@ bool KDTree::rayIntersect(const Ray &ray, Intersection &its) const {
 				const Shape *shape = m_shapes[cache->shapeIndex];
 				if (m_triangleFlag[cache->shapeIndex]) {
 					const TriMesh *trimesh = static_cast<const TriMesh *>(shape);
-					const Triangle &tri = trimesh->getTriangles()[cache->index];
+					const Triangle &tri = trimesh->getTriangles()[cache->primIndex];
 					const Point *vertexPositions = trimesh->getVertexPositions();
 					const Normal *vertexNormals = trimesh->getVertexNormals();
 					const Point2 *vertexTexcoords = trimesh->getVertexTexcoords();
@@ -195,7 +195,8 @@ bool KDTree::rayIntersect(const Ray &ray, Intersection &its) const {
 					return true;
 				} else {
 					shape->fillIntersectionRecord(ray, its.t,
-						reinterpret_cast<const uint8_t*>(temp) + 4, its);
+						reinterpret_cast<const uint8_t*>(temp) + 8, its);
+					return true;
 				}
 			}
 		}
@@ -231,7 +232,7 @@ static StatsCounter coherentPackets("General", "Coherent ray packets");
 static StatsCounter incoherentPackets("General", "Incoherent ray packets");
 
 void KDTree::rayIntersectPacket(const RayPacket4 &packet, 
-		const RayInterval4 &rayInterval, Intersection4 &its) const {
+		const RayInterval4 &rayInterval, Intersection4 &its, void *temp) const {
 	CoherentKDStackEntry MM_ALIGN16 stack[MTS_KD_MAXDEPTH];
 	RayInterval4 MM_ALIGN16 interval;
 
@@ -248,9 +249,9 @@ void KDTree::rayIntersectPacket(const RayPacket4 &packet,
 	interval.mint.ps = _mm_max_ps(interval.mint.ps, rayInterval.mint.ps);
 	interval.maxt.ps = _mm_min_ps(interval.maxt.ps, rayInterval.maxt.ps);
 
-	__m128 itsFound = _mm_cmpgt_ps(interval.mint.ps, interval.maxt.ps),
+	SSEVector itsFound( _mm_cmpgt_ps(interval.mint.ps, interval.maxt.ps));
-		   masked = itsFound;
+	__m128 masked = itsFound.ps;
-	if (_mm_movemask_ps(itsFound) == 0xF)
+	if (_mm_movemask_ps(itsFound.ps) == 0xF)
 		return;
 
 	while (currNode != NULL) {
@@ -295,54 +296,57 @@ void KDTree::rayIntersectPacket(const RayPacket4 &packet,
 		const index_type primEnd = currNode->getPrimEnd();
 
 		if (EXPECT_NOT_TAKEN(primStart != primEnd)) {
-#ifdef MTS_USE_TRIACCEL4
+			SSEVector 
-			const int count = m_packedTriangles[primStart].indirectionCount;
+				searchStart(_mm_max_ps(rayInterval.mint.ps, 
-			primStart = m_packedTriangles[primStart].indirectionIndex;
+					_mm_mul_ps(interval.mint.ps, SSEConstants::om_eps.ps))),
-			primEnd = primStart + count;
+				searchEnd(_mm_min_ps(rayInterval.maxt.ps, 
-#endif
+					_mm_mul_ps(interval.maxt.ps, SSEConstants::op_eps.ps)));
-			__m128 
-				searchStart = _mm_max_ps(rayInterval.mint.ps, 
-					_mm_mul_ps(interval.mint.ps, SSEConstants::om_eps.ps)),
-				searchEnd   = _mm_min_ps(rayInterval.maxt.ps, 
-					_mm_mul_ps(interval.maxt.ps, SSEConstants::op_eps.ps));
 
 			for (index_type entry=primStart; entry != primEnd; entry++) {
 				const TriAccel &kdTri = m_triAccel[m_indices[entry]];
 				if (EXPECT_TAKEN(kdTri.k != KNoTriangleFlag)) {
-					itsFound = _mm_or_ps(itsFound, 
+					itsFound.ps = _mm_or_ps(itsFound.ps, 
-						kdTri.rayIntersectPacket(packet, searchStart, searchEnd, masked, its));
+						kdTri.rayIntersectPacket(packet, searchStart.ps, searchEnd.ps, masked, its));
 				} else {
-#if 0
+					const Shape *shape = m_shapes[kdTri.shapeIndex];
-					/* Not a triangle - invoke the shape's intersection routine */
+
-					__m128 hasIts = m_shapes[kdTri.shapeIndex]->rayIntersectPacket(packet, 
+					for (int i=0; i<4; ++i) {
-							searchStart, searchEnd, masked, its);
+						Ray ray;
-					itsFound = _mm_or_ps(itsFound, hasIts);
+						for (int axis=0; axis<3; axis++) {
-					its.primIndex.pi  = mux_epi32(pstoepi32(hasIts), 
+							ray.o[axis] = packet.o[axis].f[i];
-						load1_epi32(kdTri.index), its.primIndex.pi);
+							ray.d[axis] = packet.d[axis].f[i];
-					its.shapeIndex.pi = mux_epi32(pstoepi32(hasIts),
+						}
-						load1_epi32(kdTri.shapeIndex), its.shapeIndex.pi);
+						Float t;
-#endif
+
+						if (shape->rayIntersect(ray, searchStart.f[i], searchEnd.f[i], t, 
+								reinterpret_cast<uint8_t *>(temp)
+								+ i * MTS_KD_INTERSECTION_TEMP + 8)) {
+							its.t.f[i] = t;
+							its.shapeIndex.i[i] = kdTri.shapeIndex;
+							its.primIndex.i[i] = KNoTriangleFlag;
+							itsFound.i[i] = 0xFFFFFFFF;
+						}
+					}
 				}
-				searchEnd = _mm_min_ps(searchEnd, its.t.ps);
+				searchEnd.ps = _mm_min_ps(searchEnd.ps, its.t.ps);
 			}
 		}
 
 		/* Abort if the tree has been traversed or if
 		   intersections have been found for all four rays */
-		if (_mm_movemask_ps(itsFound) == 0xF || --stackIndex < 0)
+		if (_mm_movemask_ps(itsFound.ps) == 0xF || --stackIndex < 0)
 			break;
 
 		/* Pop from the stack */
 		currNode = stack[stackIndex].node;
 		interval = stack[stackIndex].interval;
-		masked = _mm_or_ps(itsFound, 
+		masked = _mm_or_ps(itsFound.ps, 
 			_mm_cmpgt_ps(interval.mint.ps, interval.maxt.ps));
 	}
 }
 	
 void KDTree::rayIntersectPacketIncoherent(const RayPacket4 &packet, 
-		const RayInterval4 &rayInterval, Intersection4 &its4) const {
+		const RayInterval4 &rayInterval, Intersection4 &its4, void *temp) const {
-	uint8_t temp[MTS_KD_INTERSECTION_TEMP];
 
 	++incoherentPackets;
 	for (int i=0; i<4; i++) {
@@ -354,12 +358,13 @@ void KDTree::rayIntersectPacketIncoherent(const RayPacket4 &packet,
 		}
 		ray.mint = rayInterval.mint.f[i];
 		ray.maxt = rayInterval.maxt.f[i];
-		if (ray.mint < ray.maxt && rayIntersectHavran<false>(ray, ray.mint, ray.maxt, its4.t.f[i], temp)) {
+		if (ray.mint < ray.maxt && rayIntersectHavran<false>(ray, ray.mint, 
+				ray.maxt, its4.t.f[i], reinterpret_cast<uint8_t *>(temp) + i * MTS_KD_INTERSECTION_TEMP)) {
 			const IntersectionCache *cache = reinterpret_cast<const IntersectionCache *>(temp);
 			its4.u.f[i] = cache->u;
-			its4.v.f[i] = cache->u;
+			its4.v.f[i] = cache->v;
 			its4.shapeIndex.i[i] = cache->shapeIndex;
-			its4.primIndex.i[i] = cache->index;
+			its4.primIndex.i[i] = cache->primIndex;
 		}
 	}
 }

src/librender/preview.cpp

@@ -137,6 +137,7 @@ void PreviewWorker::processCoherent(const WorkUnit *workUnit, WorkResult *workRe
 	const SSEVector MM_ALIGN16 yOffset(0.0f, 0.0f, 1.0f, 1.0f);
 	const int pixelOffset[] = {0, 1, width, width+1};
 	const __m128 clamping = _mm_set1_ps(1/(m_minDist*m_minDist));
+	uint8_t temp[MTS_KD_INTERSECTION_TEMP*4];
 
 	const __m128 camTL[3] = {
 		 _mm_set1_ps(m_cameraTL.x),
@@ -232,9 +233,9 @@ void PreviewWorker::processCoherent(const WorkUnit *workUnit, WorkResult *workRe
 				primRay4.signs[0][0] = primSignsX ? 1 : 0;
 				primRay4.signs[1][0] = primSignsY ? 1 : 0;
 				primRay4.signs[2][0] = primSignsZ ? 1 : 0;
-				m_kdtree->rayIntersectPacket(primRay4, itv4, its4);
+				m_kdtree->rayIntersectPacket(primRay4, itv4, its4, temp);
 			} else {
-				m_kdtree->rayIntersectPacketIncoherent(primRay4, itv4, its4);
+				m_kdtree->rayIntersectPacketIncoherent(primRay4, itv4, its4, temp);
 			}
 			numRays += 4;
 
@@ -288,7 +289,6 @@ void PreviewWorker::processCoherent(const WorkUnit *workUnit, WorkResult *workRe
 				const Shape *shape = (*m_shapes)[its4.shapeIndex.i[idx]];
 				const BSDF *bsdf = shape->getBSDF();
 
-
 				if (EXPECT_TAKEN(primIndex != KNoTriangleFlag)) {
 					const TriMesh *mesh = static_cast<const TriMesh *>(shape);
 					const Triangle &t = mesh->getTriangles()[primIndex];
@@ -343,13 +343,12 @@ void PreviewWorker::processCoherent(const WorkUnit *workUnit, WorkResult *workRe
 						its.dpdv = t0.dpdv * alpha + t1.dpdv * beta + t2.dpdv * gamma;
 					}
 				} else {
-#if 0
 					Ray ray(
 						Point(primRay4.o[0].f[idx], primRay4.o[1].f[idx], primRay4.o[2].f[idx]),
 						Vector(primRay4.d[0].f[idx], primRay4.d[1].f[idx], primRay4.d[2].f[idx])
 					);
-					shape->rayIntersect(ray, its);
+					shape->fillIntersectionRecord(ray, its4.t.f[idx], temp + 
-#endif
+								+ idx * MTS_KD_INTERSECTION_TEMP + 8, its);
 				}
 
 				wo.x = nSecD[0].f[idx]; wo.y = nSecD[1].f[idx]; wo.z = nSecD[2].f[idx];
@@ -427,9 +426,9 @@ void PreviewWorker::processCoherent(const WorkUnit *workUnit, WorkResult *workRe
 				secRay4.signs[0][0] = secSignsX ? 1 : 0;
 				secRay4.signs[1][0] = secSignsY ? 1 : 0;
 				secRay4.signs[2][0] = secSignsZ ? 1 : 0;
-				m_kdtree->rayIntersectPacket(secRay4, secItv4, secIts4);
+				m_kdtree->rayIntersectPacket(secRay4, secItv4, secIts4, temp);
 			} else {
-				m_kdtree->rayIntersectPacketIncoherent(secRay4, secItv4, secIts4);
+				m_kdtree->rayIntersectPacketIncoherent(secRay4, secItv4, secIts4, temp);
 			}
 
 			for (int idx=0; idx<4; ++idx) {

src/librender/scene.cpp

@@ -239,7 +239,7 @@ void Scene::configure() {
 		Float maxExtents = std::max(extents.x, extents.y);
 		Float distance = maxExtents/(2.0f * std::tan(45 * .5f * M_PI/180));
 
-		props.setTransform("toWorld", Transform::translate(Vector(center.x, center.y, aabb.getMinimum().x - distance)));
+		props.setTransform("toWorld", Transform::translate(Vector(center.x, center.y, aabb.getMinimum().z - distance)));
 		props.setFloat("fov", 45.0f);
 
 		m_camera = static_cast<Camera *> (PluginManager::getInstance()->createObject(Camera::m_theClass, props));

src/shapes/sphere.cpp

@@ -29,136 +29,116 @@ MTS_NAMESPACE_BEGIN
  */
 
 class Sphere : public Shape {
-private:
-	Point m_center;
-	Float m_radius;
 public:
 	Sphere(const Properties &props) : Shape(props) {
-		m_radius = props.getFloat("radius", 1.0f); // Negative radius -> inside-out sphere
+		m_objectToWorld = props.getTransform("toWorld", Transform());
-		if (m_objectToWorld.hasScale()) 
+		m_center = m_objectToWorld(Point(0,0,0));
-			Log(EError, "The scale needs to be specified using the 'radius' parameter!");
+		/// non-uniform scales are not supported!
+		m_radius = m_objectToWorld(Vector(1,0,0)).length();
+		m_worldToObject = m_objectToWorld.inverse();
+		m_invSurfaceArea = 1/(4*M_PI*m_radius*m_radius);
 	}
-	
+
 	Sphere(Stream *stream, InstanceManager *manager) 
 			: Shape(stream, manager) {
+		m_objectToWorld = Transform(stream);
 		m_radius = stream->readFloat();
-		configure();
+		m_center = Point(stream);
+		m_worldToObject = m_objectToWorld.inverse();
+		m_invSurfaceArea = 1/(4*M_PI*m_radius*m_radius);
 	}
 
-	void configure() {
+	void serialize(Stream *stream, InstanceManager *manager) const {
-		Shape::configure();
+		Shape::serialize(stream, manager);
+		m_objectToWorld.serialize(stream);
+		stream->writeFloat(m_radius);
+		m_center.serialize(stream);
+	}
 
-		m_surfaceArea = 4*M_PI*m_radius*m_radius;
+	AABB getAABB() const {
-		m_invSurfaceArea = 1.0f / m_surfaceArea;
+		AABB aabb;
-		m_center = m_objectToWorld(Point(0,0,0));
 		Float absRadius = std::abs(m_radius);
-		m_aabb.min = m_center - Vector(absRadius, absRadius, absRadius);
+		aabb.min = m_center - Vector(absRadius);
-		m_aabb.max = m_center + Vector(absRadius, absRadius, absRadius);
+		aabb.max = m_center + Vector(absRadius);
-		m_bsphere.center = m_center;
+		return aabb;
-		m_bsphere.radius = absRadius;
 	}
 
-	bool rayIntersect(const Ray &ray, Float start, Float end, Float &t) const {
+	Float getSurfaceArea() const {
+		return 4*M_PI*m_radius*m_radius;
+	}
+
+	bool rayIntersect(const Ray &ray, Float mint, Float maxt, Float &t, void *tmp) const {
+		Vector ro = ray.o - m_center;
+
 		/* Transform into the local coordinate system and normalize */
-		double nearT, farT;
+		Float A = ray.d.x*ray.d.x + ray.d.y*ray.d.y + ray.d.z*ray.d.z;
-		const double ox = (double) ray.o.x - (double) m_center.x, 
+		Float B = 2 * (ray.d.x*ro.x + ray.d.y*ro.y + ray.d.z*ro.z);
-			oy = (double) ray.o.y - (double) m_center.y, 
+		Float C = ro.x*ro.x + ro.y*ro.y +
-			oz = (double) ray.o.z - (double) m_center.z;
+				ro.z*ro.z - m_radius*m_radius;
-		const double dx = ray.d.x, dy = ray.d.y, dz = ray.d.z;
-		const double A = dx*dx + dy*dy + dz*dz;
-		const double B = 2 * (dx*ox + dy*oy + dz*oz);
-		const double C = ox*ox + oy*oy + oz*oz - m_radius * m_radius;
 
-		if (!solveQuadraticDouble(A, B, C, nearT, farT))
+		Float nearT, farT;
+		if (!solveQuadratic(A, B, C, nearT, farT))
 			return false;
 
-		if (nearT > end || farT < start)
+		if (nearT > maxt || farT < mint)
 			return false;
-		if (nearT < start) {
+		if (nearT < mint) {
-			if (farT > end)
+			if (farT > maxt)
 				return false;
-			t = (Float) farT;		
+			t = farT;		
 		} else {
-			t = (Float) nearT;
+			t = nearT;
 		}
+
 		return true;
 	}
 
-	bool rayIntersect(const Ray &ray, Intersection &its) const {
+	void fillIntersectionRecord(const Ray &ray, Float t, 
-		if (!rayIntersect(ray, ray.mint, ray.maxt, its.t)) 
+			const void *temp, Intersection &its) const {
-			return false;
+		its.t = t;
-		its.p = ray(its.t);
+		its.p = ray(t);
-
+		Vector local = normalize(m_worldToObject(its.p - m_center));
-		Point local = m_worldToObject(its.p);
+		Float theta = std::acos(local.z);
-		Float absRadius = std::abs(m_radius);
-		Float theta = std::acos(local.z / absRadius);
 		Float phi = std::atan2(local.y, local.x);
+
 		if (phi < 0)
 			phi += 2*M_PI;
-		its.uv.x = theta / M_PI;
+
-		its.uv.y = phi / (2*M_PI);
+		its.uv.x = phi * (0.5 * INV_PI);
-		its.shape = this;
+		its.uv.y = theta * INV_PI;
+		its.dpdu = m_objectToWorld(Vector(-local.y, local.x, 0) * (2*M_PI));
 		Float zrad = std::sqrt(local.x*local.x + local.y*local.y);
 
-		its.geoFrame.n = Normal(normalize(its.p - m_center));
-		if (m_radius < 0)
-			its.geoFrame.n *= -1;
-
 		if (zrad > 0) {
-			Float invZRad = 1.0f / zrad;
+			Float invZRad = 1.0f / zrad,
-			Float cosPhi = local.x * invZRad;
+				  cosPhi = local.x * invZRad,
-			Float sinPhi = local.y * invZRad;
+				  sinPhi = local.y * invZRad;
-			its.dpdu = m_objectToWorld(Vector(-local.y, local.x, 0) * (2*M_PI));
 			its.dpdv = m_objectToWorld(Vector(local.z * cosPhi, local.z * sinPhi,
-					- absRadius * std::sin(theta)) * M_PI);
+					-std::sin(theta)) * M_PI);
-			its.geoFrame.s = normalize(its.dpdu - its.geoFrame.n
+    		its.geoFrame.s = normalize(its.dpdu);
-				* dot(its.geoFrame.n, its.dpdu));
+			its.geoFrame.t = normalize(its.dpdv);
-			its.geoFrame.t = cross(its.geoFrame.n, its.geoFrame.s);
+			its.geoFrame.n = normalize(cross(its.dpdv, its.dpdu));
 		} else {
 			// avoid a singularity
-			Float cosPhi = 0, sinPhi = 1;
+			const Float cosPhi = 0, sinPhi = 1;
 			its.dpdv = m_objectToWorld(Vector(local.z * cosPhi, local.z * sinPhi,
-				-absRadius*std::sin(theta))* M_PI);
+					-std::sin(theta)) * M_PI);
-			its.dpdu = cross(its.dpdv, its.geoFrame.n);
+    		its.geoFrame = Frame(normalize(its.p - m_center));
-			coordinateSystem(its.geoFrame.n, its.geoFrame.s, its.geoFrame.t);
+			if (m_radius < 0)
+				its.geoFrame.n *= -1;
 		}
 
  		its.shFrame = its.geoFrame;
  		its.wi = its.toLocal(-ray.d);
+		its.shape = this;
  		its.hasUVPartials = false;
-
-		return true;
 	}
 
-#if defined(MTS_SSE)
-	/* SSE-accelerated packet tracing is not supported for spheres at the moment */
-	__m128 rayIntersectPacket(const RayPacket4 &packet, const
-        __m128 start, __m128 end, __m128 inactive, Intersection4 &its) const {
-		SSEVector result(_mm_setzero_ps()), mint(start), maxt(end), mask(inactive);
-		Float t;
-
-		for (int i=0; i<4; i++) {
-			Ray ray;
-			for (int axis=0; axis<3; axis++) {
-				ray.o[axis] = packet.o[axis].f[i];
-				ray.d[axis] = packet.d[axis].f[i];
-			}
-			if (mask.i[i] != 0)
-				continue;
-			if (rayIntersect(ray, mint.f[i], maxt.f[i], t)) {
-				result.i[i] = 0xFFFFFFFF;
-				its.t.f[i] = t;
-			}
-		}
-		return result.ps;
-	}
-#endif
-
 	Float sampleArea(ShapeSamplingRecord &sRec, const Point2 &sample) const {
 		Vector v = squareToSphere(sample);
-		sRec.n = m_objectToWorld(Normal(v));
+		sRec.n = Normal(v);
-		sRec.p = m_objectToWorld(Point(v * m_radius));
+		sRec.p = Point(v * m_radius) + m_center;
-		return m_invSurfaceArea;
+		return 1.0f / (4*M_PI*m_radius*m_radius);
 	}
 
 	/**
@@ -193,7 +173,7 @@ public:
 
 		Ray ray(p, d);
 		Float t;
-		if (!rayIntersect(ray, 0, std::numeric_limits<Float>::infinity(), t)) {
+		if (!rayIntersect(ray, 0, std::numeric_limits<Float>::infinity(), t, NULL)) {
 			// This can happen sometimes due to roundoff errors - just fail to 
 			// generate a sample in this case.
 			return 0;
@@ -223,28 +203,25 @@ public:
 		return squareToConePdf(cosThetaMax);
 	}
 
-	void serialize(Stream *stream, InstanceManager *manager) const {
-		Shape::serialize(stream, manager);
-		stream->writeFloat(m_radius);
-	}
-
 	std::string toString() const {
 		std::ostringstream oss;
 		oss << "Sphere[" << endl
 			<< "  radius = " << m_radius << ", " << endl
 			<< "  center = " << m_center.toString() << ", " << endl
-			<< "  objectToWorld = " << indent(m_objectToWorld.toString()) << "," << endl
-			<< "  aabb = " << m_aabb.toString() << "," << endl
-			<< "  bsphere = " << m_bsphere.toString() << "," << endl
 			<< "  bsdf = " << indent(m_bsdf.toString()) << "," << endl
 			<< "  luminaire = " << indent(m_luminaire.toString()) << "," << endl
-			<< "  subsurface = " << indent(m_subsurface.toString()) << "," << endl
+			<< "  subsurface = " << indent(m_subsurface.toString()) << endl
-			<< "  surfaceArea = " << m_surfaceArea << endl
 			<< "]";
 		return oss.str();
 	}
 
 	MTS_DECLARE_CLASS()
+private:
+	Transform m_objectToWorld;
+	Transform m_worldToObject;
+	Point m_center;
+	Float m_radius;
+	Float m_invSurfaceArea;
 };
 
 MTS_IMPLEMENT_CLASS_S(Sphere, false, Shape)