diff --git a/main.go b/main.go
index a746127..6ecb5f5 100644
--- a/main.go
+++ b/main.go
@@ -6,10 +6,8 @@ import (
 	"io/ioutil"
 	"log"
 	"os"
-	"os/signal"
 	"path/filepath"
 	"runtime"
-	"runtime/pprof"
 	"strings"
 
 	. "github.com/tatsy/gopt/src/accelerator"
@@ -41,22 +39,22 @@ type JsonParam struct {
 }
 
 func main() {
-	cpuprofile := "profile.prof"
-	f, err := os.Create(cpuprofile)
-	if err != nil {
-		log.Fatal(err)
-	}
-	pprof.StartCPUProfile(f)
-	defer pprof.StopCPUProfile()
-	c := make(chan os.Signal, 1)
-	signal.Notify(c, os.Interrupt)
-	go func() {
-		for sig := range c {
-			log.Printf("captured %v, stopping profiler and exiting...", sig)
-			pprof.StopCPUProfile()
-			os.Exit(1)
-		}
-	}()
+	// cpuprofile := "profile.prof"
+	// f, err := os.Create(cpuprofile)
+	// if err != nil {
+	// 	log.Fatal(err)
+	// }
+	// pprof.StartCPUProfile(f)
+	// defer pprof.StopCPUProfile()
+	// c := make(chan os.Signal, 1)
+	// signal.Notify(c, os.Interrupt)
+	// go func() {
+	// 	for sig := range c {
+	// 		log.Printf("captured %v, stopping profiler and exiting...", sig)
+	// 		pprof.StopCPUProfile()
+	// 		os.Exit(1)
+	// 	}
+	// }()
 
 	// Parse command line args
 	var jsonFile string
diff --git a/src/bsdf/bxdf.go b/src/bsdf/bxdf.go
new file mode 100644
index 0000000..0256fa9
--- /dev/null
+++ b/src/bsdf/bxdf.go
@@ -0,0 +1,5 @@
+package bsdf
+
+const (
+	deltaEps = 1.0e-12
+)
diff --git a/src/bsdf/lambert.go b/src/bsdf/lambert.go
index 2ec088d..d5b5e1c 100644
--- a/src/bsdf/lambert.go
+++ b/src/bsdf/lambert.go
@@ -21,7 +21,10 @@ func (f *LambertReflection) Eval(wi, wo *Vector3d) *Color {
 }
 
 func (f *LambertReflection) Pdf(wi, wo *Vector3d) Float {
-	return 1.0
+	if SameHemisphere(wi, wo) {
+		return AbsCosTheta(wi) / math.Pi
+	}
+	return 0.0
 }
 
 func (f *LambertReflection) Sample(wo *Vector3d, u *Vector2d) (*Color, *Vector3d, Float) {
diff --git a/src/bsdf/specular.go b/src/bsdf/specular.go
index 0507683..a7c4a90 100644
--- a/src/bsdf/specular.go
+++ b/src/bsdf/specular.go
@@ -16,14 +16,14 @@ func NewSpecularReflection(re *Color) *SpecularReflection {
 }
 
 func (f *SpecularReflection) Eval(wi, wo *Vector3d) *Color {
-	if NewVector3d(-wi.X, -wi.Y, wi.Z).Dot(wo) >= 1.0-Eps {
-		return f.re
+	if NewVector3d(-wi.X, -wi.Y, wi.Z).Dot(wo) >= 1.0-deltaEps {
+		return f.re.Scale(1.0 / AbsCosTheta(wi))
 	}
 	return NewColor(0.0, 0.0, 0.0)
 }
 
 func (f *SpecularReflection) Pdf(wi, wo *Vector3d) Float {
-	if NewVector3d(-wi.X, -wi.Y, wi.Z).Dot(wo) >= 1.0-Eps {
+	if NewVector3d(-wi.X, -wi.Y, wi.Z).Dot(wo) >= 1.0-deltaEps {
 		return 1.0
 	}
 	return 0.0
@@ -58,7 +58,7 @@ func (f *SpecularFresnel) Eval(wi, wo *Vector3d) *Color {
 	cosThetaI := CosTheta(wo)
 	F := FresnelDielectric(cosThetaI, f.etaA, f.etaB)
 	if SameHemisphere(wi, wo) {
-		if NewVector3d(-wi.X, -wi.Y, wi.Z).Dot(wo) >= 1.0-Eps {
+		if NewVector3d(-wi.X, -wi.Y, wi.Z).Dot(wo) >= 1.0-deltaEps {
 			return f.re.Scale(F / AbsCosTheta(wi))
 		}
 	} else {
@@ -74,7 +74,7 @@ func (f *SpecularFresnel) Eval(wi, wo *Vector3d) *Color {
 		}
 
 		wt, isRefr := Refract(wo, faceForwardNormal, etaI/etaT)
-		if isRefr && wt.Dot(wi) > 1.0-Eps {
+		if isRefr && wt.Dot(wi) > 1.0-deltaEps {
 			return f.tr.Scale((1.0 - F) / AbsCosTheta(wi))
 		}
 	}
@@ -85,7 +85,7 @@ func (f *SpecularFresnel) Pdf(wi, wo *Vector3d) Float {
 	cosThetaI := CosTheta(wo)
 	F := FresnelDielectric(cosThetaI, f.etaA, f.etaB)
 	if SameHemisphere(wi, wo) {
-		if NewVector3d(-wi.X, -wi.Y, wi.Z).Dot(wo) >= 1.0-Eps {
+		if NewVector3d(-wi.X, -wi.Y, wi.Z).Dot(wo) >= 1.0-deltaEps {
 			return F
 		}
 	} else {
@@ -101,7 +101,7 @@ func (f *SpecularFresnel) Pdf(wi, wo *Vector3d) Float {
 		}
 
 		wt, isRefract := Refract(wo, faceForwardNormal, etaI/etaT)
-		if isRefract && wt.Dot(wi) > 1.0-Eps {
+		if isRefract && wt.Dot(wi) > 1.0-deltaEps {
 			return 1.0 - F
 		}
 	}
diff --git a/src/core/bsdf.go b/src/core/bsdf.go
index 3b1d1ac..7acae9f 100644
--- a/src/core/bsdf.go
+++ b/src/core/bsdf.go
@@ -32,7 +32,19 @@ func (bsdf *Bsdf) Eval(wi, wo *Vector3d) *Color {
 	zLocal = bsdf.ns.Dot(wo)
 	woLocal := NewVector3d(xLocal, yLocal, zLocal)
 	return bsdf.bxdf.Eval(wiLocal, woLocal)
+}
 
+func (bsdf *Bsdf) Pdf(wi, wo *Vector3d) Float {
+	var xLocal, yLocal, zLocal Float
+	xLocal = bsdf.ts.Dot(wi)
+	yLocal = bsdf.bs.Dot(wi)
+	zLocal = bsdf.ns.Dot(wi)
+	wiLocal := NewVector3d(xLocal, yLocal, zLocal)
+	xLocal = bsdf.ts.Dot(wo)
+	yLocal = bsdf.bs.Dot(wo)
+	zLocal = bsdf.ns.Dot(wo)
+	woLocal := NewVector3d(xLocal, yLocal, zLocal)
+	return bsdf.bxdf.Pdf(wiLocal, woLocal)
 }
 
 func (bsdf *Bsdf) SampleWi(wo *Vector3d, u *Vector2d) (*Color, *Vector3d, Float, int) {
diff --git a/src/core/color.go b/src/core/color.go
index 7343c84..062eaad 100644
--- a/src/core/color.go
+++ b/src/core/color.go
@@ -1,59 +1,71 @@
 package core
 
 import (
-    "fmt"
+	"fmt"
+	"math"
 )
 
 type Color struct {
-    R, G, B Float
+	R, G, B Float
 }
 
 func NewColor(r, g, b Float) *Color {
-    c := &Color{}
-    c.R = r
-    c.G = g
-    c.B = b
-    return c
+	c := &Color{}
+	c.R = r
+	c.G = g
+	c.B = b
+	return c
 }
 
 func NewColorWithString(s string) *Color {
-    var r, g, b Float
-    n, _ := fmt.Sscanf(s, "(%f, %f, %f)", &r, &g, &b)
-    if n != 3 {
-        panic(fmt.Sprintf("Failed to parse Vector3d: %s", s))
-    }
-    return &Color{r, g, b}
+	var r, g, b Float
+	n, _ := fmt.Sscanf(s, "(%f, %f, %f)", &r, &g, &b)
+	if n != 3 {
+		panic(fmt.Sprintf("Failed to parse Vector3d: %s", s))
+	}
+	return &Color{r, g, b}
 }
 
-
 func (c *Color) Y() Float {
-    return 0.299 * c.R + 0.587 * c.G + 0.114 * c.B
+	return 0.299*c.R + 0.587*c.G + 0.114*c.B
 }
 
 func (c *Color) IsBlack() bool {
-    return c.R == 0.0 && c.G == 0.0 && c.B == 0.0
+	return c.R == 0.0 && c.G == 0.0 && c.B == 0.0
+}
+
+func (c *Color) IsNaN() bool {
+	return math.IsNaN(c.R) || math.IsNaN(c.G) || math.IsNaN(c.B)
+}
+
+func (c *Color) IsInf() bool {
+	return math.IsInf(c.R, 0) || math.IsInf(c.G, 0) || math.IsInf(c.B, 0)
+}
+
+func (c *Color) IsValid() bool {
+	return !c.IsInf() && !c.IsNaN()
 }
 
 func (c1 *Color) Add(c2 *Color) *Color {
-    ret := &Color{}
-    ret.R = c1.R + c2.R
-    ret.G = c1.G + c2.G
-    ret.B = c1.B + c2.B
-    return ret
+	ret := &Color{}
+	ret.R = c1.R + c2.R
+	ret.G = c1.G + c2.G
+	ret.B = c1.B + c2.B
+	return ret
 }
 
 func (c *Color) Scale(x Float) *Color {
-    ret := &Color{}
-    ret.R = c.R * x
-    ret.G = c.G * x
-    ret.B = c.B * x
-    return ret
+	ret := &Color{}
+	ret.R = c.R * x
+	ret.G = c.G * x
+	ret.B = c.B * x
+	return ret
 }
 
 func (c1 *Color) Multiply(c2 *Color) *Color {
-    ret := &Color{}
-    ret.R = c1.R * c2.R
-    ret.G = c1.G * c2.G
-    ret.B = c1.B * c2.B
-    return ret
+	ret := &Color{}
+	ret.R = c1.R * c2.R
+	ret.G = c1.G * c2.G
+	ret.B = c1.B * c2.B
+	return ret
 }
diff --git a/src/core/interface.go b/src/core/interface.go
index d738b5f..04bbc22 100644
--- a/src/core/interface.go
+++ b/src/core/interface.go
@@ -6,7 +6,8 @@ type Accelerator interface {
 
 type Shape interface {
 	Intersect(ray *Ray, tHit *Float, isect *Intersection) bool
-	SampleP(u *Vector2d) (*Vector3d, *Vector3d, Float)
+	SamplePoint(u *Vector2d) (*Vector3d, *Vector3d)
+	Pdf(isect *Intersection, wi *Vector3d) Float
 	Bounds() *Bounds3d
 }
 
@@ -32,4 +33,5 @@ type Light interface {
 	Le() *Color
 	LeWithRay(ray *Ray) *Color
 	SampleLi(isect *Intersection, u *Vector2d) (*Color, *Vector3d, Float, *VisibilityTester)
+	PdfLi(isect *Intersection, wi *Vector3d) Float
 }
diff --git a/src/core/intersection.go b/src/core/intersection.go
index 719b101..a8a21e9 100644
--- a/src/core/intersection.go
+++ b/src/core/intersection.go
@@ -1,31 +1,33 @@
 package core
 
+import "math"
+
 type Intersection struct {
-    Pos, Normal, Wo *Vector3d
-    primitive *Primitive
+	Pos, Normal, Wo *Vector3d
+	primitive       *Primitive
 }
 
 func NewIntersection(pos, normal, wo *Vector3d) *Intersection {
-    isect := &Intersection{}
-    isect.Pos = pos
-    isect.Normal = normal
-    isect.Wo = wo
-    isect.primitive = nil
-    return isect
+	isect := &Intersection{}
+	isect.Pos = pos
+	isect.Normal = normal
+	isect.Wo = wo
+	isect.primitive = nil
+	return isect
 }
 
 func (isect *Intersection) SpawnRay(wi *Vector3d) *Ray {
-    return NewRay(isect.Pos, wi)
+	return NewRay(isect.Pos, wi)
 }
 
 func (isect *Intersection) Bsdf() *Bsdf {
-    return NewBsdf(isect, isect.primitive.Bxdf())
+	return NewBsdf(isect, isect.primitive.Bxdf())
 }
 
 func (isect *Intersection) Le(wi *Vector3d) *Color {
-    if isect.primitive.Light() == nil {
-        return NewColor(0.0, 0.0, 0.0)
-    }
-    dot := wi.Dot(isect.Normal)
-    return isect.primitive.Light().Le().Scale(dot)
+	if isect.primitive.Light() == nil {
+		return NewColor(0.0, 0.0, 0.0)
+	}
+	dot := math.Max(0.0, wi.Dot(isect.Normal))
+	return isect.primitive.Light().Le().Scale(dot)
 }
diff --git a/src/core/ray.go b/src/core/ray.go
index 23d4356..8610168 100644
--- a/src/core/ray.go
+++ b/src/core/ray.go
@@ -1,65 +1,66 @@
 package core
 
 import (
-    "math"
+	"math"
 )
 
+// Ray is exhibited from the position Org to the direction Dir.
 type Ray struct {
-    Org *Vector3d
-    Dir *Vector3d
-    InvDir *Vector3d
-    MaxDist Float
+	Org     *Vector3d
+	Dir     *Vector3d
+	InvDir  *Vector3d
+	MaxDist Float
 }
 
 func NewRay(org *Vector3d, dir *Vector3d) *Ray {
-    r := &Ray{}
-    r.Org = org
-    r.Dir = dir
-    r.MaxDist = Infinity
-    r.InvDir = InvertDir(r.Dir)
-    return r
+	r := new(Ray)
+	r.Org = org
+	r.Dir = dir
+	r.MaxDist = Infinity
+	r.InvDir = invertDir(r.Dir)
+	return r
 }
 
 func NewRayBetweenPoints(origin, target *Vector3d) *Ray {
-    dir := target.Subtract(origin)
-    dist := dir.Length()
-    dir = dir.Divide(dist)
+	dir := target.Subtract(origin)
+	dist := dir.Length()
+	dir = dir.Divide(dist)
 
-    ray := &Ray{}
-    ray.Org = origin.Add(dir.Scale(Eps))
-    ray.Dir = dir
-    ray.InvDir = InvertDir(dir)
-    ray.MaxDist = dist - 2.0 * Eps
-    return ray
+	ray := &Ray{}
+	ray.Org = origin.Add(dir.Scale(Eps))
+	ray.Dir = dir
+	ray.InvDir = invertDir(dir)
+	ray.MaxDist = dist - 2.0*Eps
+	return ray
 }
 
 func (ray *Ray) Clone() *Ray {
-    return &Ray{
-        Org: ray.Org,
-        Dir: ray.Dir,
-        InvDir: ray.InvDir,
-        MaxDist: ray.MaxDist,
-    }
+	return &Ray{
+		Org:     ray.Org,
+		Dir:     ray.Dir,
+		InvDir:  ray.InvDir,
+		MaxDist: ray.MaxDist,
+	}
 }
 
-func InvertDir(v *Vector3d) *Vector3d {
-    d := &Vector3d{}
-    if (math.Abs(v.X) > Eps) {
-        d.X = 1.0 / v.X
-    } else {
-        d.X = Infinity * Sign(v.X)
-    }
+func invertDir(v *Vector3d) *Vector3d {
+	d := &Vector3d{}
+	if math.Abs(v.X) > Eps {
+		d.X = 1.0 / v.X
+	} else {
+		d.X = Infinity * Sign(v.X)
+	}
 
-    if (math.Abs(v.Y) > Eps) {
-        d.Y = 1.0 / v.Y
-    } else {
-        d.Y = Infinity * Sign(v.Y)
-    }
+	if math.Abs(v.Y) > Eps {
+		d.Y = 1.0 / v.Y
+	} else {
+		d.Y = Infinity * Sign(v.Y)
+	}
 
-    if (math.Abs(v.Z) > Eps) {
-        d.Z = 1.0 / v.Z
-    } else {
-        d.Z = Infinity * Sign(v.Z)
-    }
-    return d
+	if math.Abs(v.Z) > Eps {
+		d.Z = 1.0 / v.Z
+	} else {
+		d.Z = Infinity * Sign(v.Z)
+	}
+	return d
 }
diff --git a/src/core/vector3d.go b/src/core/vector3d.go
index 3a64380..7155387 100644
--- a/src/core/vector3d.go
+++ b/src/core/vector3d.go
@@ -5,26 +5,32 @@ import (
 	"math"
 )
 
+// Vector3d is a 3D vector.
 type Vector3d struct {
 	X, Y, Z Float
 }
 
+// NewVector3d creates a new vector with specified coordinates.
 func NewVector3d(x, y, z Float) *Vector3d {
-	ret := &Vector3d{x, y, z}
+	ret := new(Vector3d)
+	ret.X = x
+	ret.Y = y
+	ret.Z = z
 	return ret
 }
 
+// NewVector3dWithString parses a string, and returns a parsed vector.
 func NewVector3dWithString(s string) *Vector3d {
 	var x, y, z Float
 	n, _ := fmt.Sscanf(s, "(%f, %f, %f)", &x, &y, &z)
 	if n != 3 {
 		panic(fmt.Sprintf("Failed to parse Vector3d: %s", s))
 	}
-	return &Vector3d{x, y, z}
+	return NewVector3d(x, y, z)
 }
 
 func (v1 *Vector3d) Add(v2 *Vector3d) *Vector3d {
-	ret := &Vector3d{}
+	ret := new(Vector3d)
 	ret.X = v1.X + v2.X
 	ret.Y = v1.Y + v2.Y
 	ret.Z = v1.Z + v2.Z
@@ -32,7 +38,7 @@ func (v1 *Vector3d) Add(v2 *Vector3d) *Vector3d {
 }
 
 func (v *Vector3d) Negate() *Vector3d {
-	ret := &Vector3d{}
+	ret := new(Vector3d)
 	ret.X = -v.X
 	ret.Y = -v.Y
 	ret.Z = -v.Z
@@ -44,7 +50,7 @@ func (v1 *Vector3d) Subtract(v2 *Vector3d) *Vector3d {
 }
 
 func (v1 *Vector3d) Scale(s Float) *Vector3d {
-	ret := &Vector3d{}
+	ret := new(Vector3d)
 	ret.X = v1.X * s
 	ret.Y = v1.Y * s
 	ret.Z = v1.Z * s
diff --git a/src/core/vector3d_test.go b/src/core/vector3d_test.go
index cba7b39..8e02276 100644
--- a/src/core/vector3d_test.go
+++ b/src/core/vector3d_test.go
@@ -32,16 +32,26 @@ func TestVector3dInistance(t *testing.T) {
 }
 
 func TestVector3dConstructWithString(t *testing.T) {
-	testCases := map[string][3]Float{
+	testCases := map[string][]Float{
 		"(1, 2, 3)":    {1.0, 2.0, 3.0},
 		"(0, 0, 0)":    {0.0, 0.0, 0.0},
 		"(-1, -2, -3)": {-1.0, -2.0, -3.0},
+		"(-1, -2, A)":  {},
 	}
 
 	for s, vec := range testCases {
-		v := NewVector3dWithString(s)
-		if v.X != vec[0] || v.Y != vec[1] || v.Z != vec[2] {
-			t.Errorf("%s converted to %v", s, v)
+		if len(vec) != 0 {
+			v := NewVector3dWithString(s)
+			if v.X != vec[0] || v.Y != vec[1] || v.Z != vec[2] {
+				t.Errorf("%s converted to %v", s, v)
+			}
+		} else {
+			defer func() {
+				if p := recover(); p != nil {
+				}
+			}()
+			NewVector3dWithString(s)
+			t.Errorf("String \"%s\" must not be parsed.", s)
 		}
 	}
 }
diff --git a/src/integrator/path.go b/src/integrator/path.go
index 9a5e04f..37b26c2 100644
--- a/src/integrator/path.go
+++ b/src/integrator/path.go
@@ -105,19 +105,55 @@ func NextEventEstimation(isect *Intersection, scene *Scene, sampler Sampler) *Co
 		return NewColor(0.0, 0.0, 0.0)
 	}
 
-	numLights := len(scene.Lights)
-	lightId := int(sampler.Get1D() * Float(numLights))
-	light := scene.Lights[lightId]
+	Ld := NewColor(0.0, 0.0, 0.0)
+	lightPdf, bsdfPdf := 0.0, 0.0
 
-	Li, wi, pdf, vis := light.SampleLi(isect, sampler.Get2D())
-	if pdf == 0.0 {
-		return NewColor(0.0, 0.0, 0.0)
+	// Light sampling
+	numLights := len(scene.Lights)
+	lightID := int(sampler.Get1D() * Float(numLights))
+	light := scene.Lights[lightID]
+
+	Li, wi, lightPdf, vis := light.SampleLi(isect, sampler.Get2D())
+	if !Li.IsBlack() && lightPdf > 0.0 && vis.Unoccluded(scene) {
+		bsdfPdf = isect.Bsdf().Pdf(wi, isect.Wo)
+		weight := lightPdf / (lightPdf + bsdfPdf)
+		f := isect.Bsdf().Eval(wi, isect.Wo).Scale(math.Abs(isect.Normal.Dot(wi)))
+		L := f.Multiply(Li).Scale(weight / lightPdf)
+		Ld = Ld.Add(L)
 	}
 
-	if !vis.Unoccluded(scene) {
-		return NewColor(0.0, 0.0, 0.0)
+	// BSDF sampling
+	lightPdf, bsdfPdf = 0.0, 0.0
+	f, wi, bsdfPdf, sampledType := isect.Bsdf().SampleWi(isect.Wo, sampler.Get2D())
+	if !f.IsBlack() && bsdfPdf > 0.0 {
+		sampledSpecular := (sampledType & BSDF_SPECULAR) != 0
+		f = f.Scale(math.Abs(isect.Normal.Dot(wi)))
+		weight := 1.0
+		if !sampledSpecular {
+			lightPdf = light.PdfLi(isect, wi)
+			if lightPdf <= 0.0 {
+				return Ld.Scale(Float(numLights))
+			}
+			weight = bsdfPdf / (bsdfPdf + lightPdf)
+		}
+
+		var testIsect Intersection
+		ray := NewRay(isect.Pos, wi)
+		isIntersect := scene.Intersect(ray, &testIsect)
+		Li := NewColor(0.0, 0.0, 0.0)
+		if isIntersect {
+			rr := testIsect.Le(ray.Dir.Negate())
+			Li = Li.Add(rr)
+		} else {
+			for _, l := range scene.Lights {
+				rr := l.LeWithRay(ray)
+				Li = Li.Add(rr)
+			}
+		}
+
+		L := f.Multiply(Li).Scale(weight / bsdfPdf)
+		Ld = Ld.Add(L)
 	}
 
-	f := isect.Bsdf().Eval(wi, isect.Wo).Scale(math.Abs(isect.Normal.Dot(wi)))
-	return f.Multiply(Li).Scale(1.0 / pdf).Scale(Float(numLights))
+	return Ld.Scale(Float(numLights))
 }
diff --git a/src/light/area.go b/src/light/area.go
index 5278b81..93cdf67 100644
--- a/src/light/area.go
+++ b/src/light/area.go
@@ -20,20 +20,24 @@ func NewAreaLight(shape Shape, radiance *Color) *AreaLight {
 }
 
 func (l *AreaLight) SampleLi(isect *Intersection, u *Vector2d) (*Color, *Vector3d, Float, *VisibilityTester) {
-	pos, normal, pdf := l.shape.SampleP(u)
+	pos, normal := l.shape.SamplePoint(u)
 	vt := NewVisibilityTester(isect.Pos, pos)
 
-	wi := pos.Subtract(isect.Pos)
-	distSquared := wi.LengthSquared()
+	wi := pos.Subtract(isect.Pos).Normalized()
+	pdf := l.shape.Pdf(isect, wi)
 
-	wi = wi.Normalized()
-	dot0 := math.Max(0.0, wi.Dot(isect.Normal))
-	dot1 := math.Max(0.0, -wi.Dot(normal))
-
-	Le := l.radiance.Scale(dot0 * dot1 / distSquared)
+	dot := math.Max(0.0, -wi.Dot(normal))
+	Le := l.radiance
+	if dot <= 0.0 {
+		Le = NewColor(0.0, 0.0, 0.0)
+	}
 	return Le, wi, pdf, vt
 }
 
+func (l *AreaLight) PdfLi(isect *Intersection, wi *Vector3d) Float {
+	return l.shape.Pdf(isect, wi)
+}
+
 func (l *AreaLight) Le() *Color {
 	return l.radiance
 }
diff --git a/src/shape/triangle.go b/src/shape/triangle.go
index 95080bc..2750c05 100644
--- a/src/shape/triangle.go
+++ b/src/shape/triangle.go
@@ -6,12 +6,15 @@ import (
 	. "github.com/tatsy/gopt/src/core"
 )
 
+// Triangle is a triangle which holds three positions,
+// normals and texture cooridnates.
 type Triangle struct {
 	Points    [3]*Vector3d
 	Normals   [3]*Vector3d
 	TexCoords [3]*Vector2d
 }
 
+// NewTriangleWithP create a triangle with three positions.
 func NewTriangleWithP(points [3]*Vector3d) *Triangle {
 	t := &Triangle{}
 	t.Points = points
@@ -99,7 +102,7 @@ func (t *Triangle) Intersect(ray *Ray, tHit *Float, isect *Intersection) bool {
 	return true
 }
 
-func (t *Triangle) SampleP(rnd *Vector2d) (*Vector3d, *Vector3d, Float) {
+func (t *Triangle) SamplePoint(rnd *Vector2d) (*Vector3d, *Vector3d) {
 	u, v := rnd.X, rnd.Y
 	if u+v >= 1.0 {
 		u = 1.0 - u
@@ -112,12 +115,27 @@ func (t *Triangle) SampleP(rnd *Vector2d) (*Vector3d, *Vector3d, Float) {
 	normal := t.Normals[0].Scale(1.0 - u - v).
 		Add(t.Normals[1].Scale(u)).
 		Add(t.Normals[2].Scale(v))
-	area := 0.5 * (t.Points[1].Subtract(t.Points[0])).Cross(t.Points[2].Subtract(t.Points[0])).Length()
-	pdf := 0.0
-	if area != 0.0 {
-		pdf = 1.0 / area
+
+	return pos, normal
+}
+
+func (t *Triangle) Pdf(isect *Intersection, wi *Vector3d) Float {
+	ray := NewRay(isect.Pos, wi)
+	var tHit Float
+	var temp Intersection
+	if !t.Intersect(ray, &tHit, &temp) {
+		return 0.0
 	}
-	return pos, normal, pdf
+
+	dot := math.Max(0.0, -temp.Normal.Dot(ray.Dir))
+	dist2 := temp.Pos.Subtract(isect.Pos).LengthSquared()
+	return dist2 / (dot * t.Area())
+}
+
+func (t *Triangle) Area() Float {
+	e1 := t.Points[1].Subtract(t.Points[0])
+	e2 := t.Points[2].Subtract(t.Points[0])
+	return 0.5 * e1.Cross(e2).Length()
 }
 
 func (t *Triangle) Bounds() *Bounds3d {