From 7b981cdb60d0ef80a688848bd5c0c661af826e46 Mon Sep 17 00:00:00 2001 From: tatsy Date: Sat, 11 Nov 2017 18:26:40 +0900 Subject: [PATCH] SAH implemented. --- Makefile | 5 +- README.md | 18 ++++ data/scene.json | 6 +- src/accelerator/bvh.go | 201 +++++++++++++++++------------------ src/accelerator/bvh_test.go | 90 ++++++++++++---- src/accelerator/bvh_utils.go | 162 ++++++++++++++++++++++++++++ src/core/bounds3d.go | 22 +++- src/core/bounds3d_test.go | 29 +++++ src/main.go | 38 ++++--- 9 files changed, 421 insertions(+), 150 deletions(-) create mode 100644 src/accelerator/bvh_utils.go diff --git a/Makefile b/Makefile index 66cf9c2..d6ddcd9 100644 --- a/Makefile +++ b/Makefile @@ -1,14 +1,13 @@ export GOPATH = $(PWD) - all: go build ./... run: - go build ./src/main.go && ./main + go build ./src/main.go && ./main ${ARGS} test: - go test ./... --cover + go test ./... --cover -v clean: rm -rf main *.jpg *.png *.prof diff --git a/README.md b/README.md index d37b452..3e8f6c4 100644 --- a/README.md +++ b/README.md @@ -3,3 +3,21 @@ gopt [![Build Status](https://travis-ci.org/tatsy/gopt.svg?branch=master)](https://travis-ci.org/tatsy/gopt) [![Coverage Status](https://coveralls.io/repos/github/tatsy/gopt/badge.svg)](https://coveralls.io/github/tatsy/gopt) + +> Physically based path tracer implemented with Go. + +## Usage + +```sh +export GOPATH=`pwd` +go build ./src/main.go +./main -i [Scene JSON file] +``` + +## Result + +gopher.jpg + +## Copyright + +MIT License 2017 (c) Tatsuya Yatagawa (tatsy) diff --git a/data/scene.json b/data/scene.json index 5c79c9a..9172a8d 100644 --- a/data/scene.json +++ b/data/scene.json @@ -2,8 +2,8 @@ { "name": "image", "parameters": [ - {"name": "width", "value": "1280"}, - {"name": "height", "value": "720"} + {"name": "width", "value": "640"}, + {"name": "height", "value": "360"} ] }, { @@ -11,7 +11,7 @@ "parameters": [ {"name": "type", "value": "path"}, {"name": "max-bounces", "value": "16"}, - {"name": "num-samples", "value": "128"} + {"name": "num-samples", "value": "32"} ] }, { diff --git a/src/accelerator/bvh.go b/src/accelerator/bvh.go index 5b387c0..4e489da 100644 --- a/src/accelerator/bvh.go +++ b/src/accelerator/bvh.go @@ -1,88 +1,10 @@ package accelerator import ( - "sort" + "math" . "core" ) -type BvhNode struct { - left, right *BvhNode - primId int - bbox *Bounds3d -} - -func NewLeafNode(primId int, b *Bounds3d) *BvhNode { - node := &BvhNode{} - node.left = nil - node.right = nil - node.primId = primId - node.bbox = b - return node -} - -func NewForkNode(left *BvhNode, right *BvhNode, b *Bounds3d) *BvhNode { - node := &BvhNode{} - node.left = left - node.right = right - node.primId = -1 - node.bbox = b - return node -} - -func (node *BvhNode) IsLeaf() bool { - return node.left == nil && node.right == nil -} - -type SortItem struct { - v *Vector3d - i int -} - -func NewSortItem(v *Vector3d, i int) *SortItem { - item := &SortItem{} - item.v = v - item.i = i - return item -} - -type AxisSorter struct { - Items []*SortItem - Axis int -} - -func NewAxisSorter(items []*SortItem, axis int) *AxisSorter { - sorter := &AxisSorter{} - sorter.Items = items - sorter.Axis = axis - return sorter -} - -func (a *AxisSorter) Len() int { - return len(a.Items) -} - -func (a *AxisSorter) Swap(i, j int) { - a.Items[i], a.Items[j] = a.Items[j], a.Items[i] -} - -func (a *AxisSorter) Less(i, j int) bool { - v1 := a.Items[i].v - v2 := a.Items[j].v - return v1.NthElement(a.Axis) < v2.NthElement(a.Axis) -} - -type IndexedPrimitive struct { - p *Primitive - i int -} - -func NewIndexedPrimitive(p *Primitive, i int) *IndexedPrimitive { - ip := &IndexedPrimitive{} - ip.p = p - ip.i = i - return ip -} - type Bvh struct { primitives []*Primitive nodes []*BvhNode @@ -93,44 +15,119 @@ func NewBvh(primitives []*Primitive) *Bvh { bvh := &Bvh{} bvh.primitives = primitives - ips := make([]*IndexedPrimitive, len(primitives)) + buildData := make([]*BvhPrimitiveInfo, len(primitives)) for i, p := range primitives { - ips[i] = NewIndexedPrimitive(p, i) + buildData[i] = NewBvhPrimitiveInfo(i, p.Bounds()) } - bvh.root = NewBvhSub(bvh, ips) + bvh.root = NewBvhSub(bvh, buildData) return bvh } -func NewBvhSub(bvh *Bvh, primitives []*IndexedPrimitive) *BvhNode { - if len(primitives) == 1 { - node := NewLeafNode(primitives[0].i, primitives[0].p.Bounds()) +func NewBvhSub(bvh *Bvh, buildData []*BvhPrimitiveInfo) *BvhNode { + numData := len(buildData) + if numData == 1 { + node := NewLeafNode(buildData[0].primitiveId, buildData[0].bounds) bvh.nodes = append(bvh.nodes, node) return node } - bbox := NewBounds3d() - items := make([]*SortItem, len(primitives)) - for i := range primitives { - b := primitives[i].p.Bounds() - bbox.Merge(b) - items[i] = &SortItem{b.Center(), i} + bounds := NewBounds3d() + for i := 0; i < numData; i++ { + bounds = bounds.Merge(buildData[i].bounds) } - axis := bbox.MaxExtent() + splitAxis := bounds.MaxExtent() + + splitMethod := BVH_SPLIT_METHOD_SAH + splitPos := numData / 2 + + switch splitMethod { + case BVH_SPLIT_METHOD_EQUAL_COUNT: + sortable := NewBvhPrimitiveSortable( + buildData, splitAxis, nil, + ) + SliceNthElement(sortable, 0, numData, splitPos) + case BVH_SPLIT_METHOD_SAH: + if numData <= 4 { + sortable := NewBvhPrimitiveSortable( + buildData, splitAxis, nil, + ) + SliceNthElement(sortable, 0, numData, splitPos) + } else { + numBuckets := 12 + centroidBounds := NewBounds3d() + for i := 0; i < numData; i++ { + centroidBounds = centroidBounds.MergePoint(buildData[i].centroid) + } - axisSorter := &AxisSorter{items, axis} - sort.Sort(axisSorter) + buckets := make([]*BucketInfo, numBuckets) + for i := 0; i < numBuckets; i++ { + buckets[i] = NewBucketInfo(0, NewBounds3d()) + } - newPrimitives := make([]*IndexedPrimitive, len(primitives)) - for i := range items { - newPrimitives[i] = primitives[items[i].i] - } + cMin := centroidBounds.MinPos.NthElement(splitAxis) + cMax := centroidBounds.MaxPos.NthElement(splitAxis) + invDenom := 1.0 / (math.Abs(cMax - cMin) + Eps) + for i := 0; i < numData; i++ { + c0 := buildData[i].centroid.NthElement(splitAxis) + c1 := centroidBounds.MinPos.NthElement(splitAxis) + numer := c0 - c1 + b := int(Float(numBuckets) * math.Abs(numer) * invDenom) + if b == numBuckets { + b = numBuckets - 1 + } + + buckets[b].count += 1 + buckets[b].bounds = buckets[b].bounds.Merge(buildData[i].bounds) + } - iHalf := len(newPrimitives) / 2 - leftNode := NewBvhSub(bvh, newPrimitives[:iHalf]) - rightNode := NewBvhSub(bvh, newPrimitives[iHalf:]) + bucketCost := make([]Float, numBuckets - 1) + for i := 0; i < numBuckets - 1; i++ { + b0 := NewBounds3d() + b1 := NewBounds3d() + cnt0, cnt1 := 0, 0 + for j := 0; j <= i; j++ { + b0 = b0.Merge(buckets[j].bounds) + cnt0 += buckets[j].count + } + for j := i + 1; j < numBuckets; j++ { + b1 = b1.Merge(buckets[j].bounds) + cnt1 += buckets[j].count + } + bucketCost[i] += 0.125 + (Float(cnt0) * b0.Area() + Float(cnt1) * b1.Area()) / bounds.Area() + } + + minCost := bucketCost[0] + minCostSplit := 0 + for i := 1; i < numBuckets - 1; i++ { + if minCost > bucketCost[i] { + minCost = bucketCost[i] + minCostSplit = i + } + } + + if minCost < Float(numData) { + comparator := func(info *BvhPrimitiveInfo) bool { + cMin := centroidBounds.MinPos.NthElement(splitAxis) + cMax := centroidBounds.MaxPos.NthElement(splitAxis) + inv := 1.0 / (math.Abs(cMax - cMin) + Eps) + diff := math.Abs(info.centroid.NthElement(splitAxis) - cMin) + b := int(Float(numBuckets) * diff * inv) + if b >= numBuckets { + b = numBuckets - 1 + } + return b <= minCostSplit + } + splitPos = SlicePartition(NewBvhPrimitiveSortable( + buildData, splitAxis, comparator, + )) + } + } + } - node := NewForkNode(leftNode, rightNode, bbox) + leftNode := NewBvhSub(bvh, buildData[:splitPos]) + rightNode := NewBvhSub(bvh, buildData[splitPos:]) + node := NewForkNode(leftNode, rightNode, bounds) bvh.nodes = append(bvh.nodes, node) return node } diff --git a/src/accelerator/bvh_test.go b/src/accelerator/bvh_test.go index 6be4794..48f3756 100644 --- a/src/accelerator/bvh_test.go +++ b/src/accelerator/bvh_test.go @@ -1,7 +1,6 @@ package accelerator import ( - "sort" "math" "math/rand" "testing" @@ -9,25 +8,80 @@ import ( . "shape" ) -func TestAxisSorter(t *testing.T) { - items := make([]*SortItem, 100) - for i := range items { - v := NewVector3d( - rand.Float64(), - rand.Float64(), - rand.Float64(), - ) - items[i] = &SortItem{v, i} +type BvhTestSotable struct { + values []int +} + +func NewBvhTestSorable(values []int) *BvhTestSotable { + s := &BvhTestSotable{} + s.values = values + return s +} + +func (s *BvhTestSotable) Len() int { + return len(s.values) +} + +func (s *BvhTestSotable) Swap(i, j int) { + s.values[i], s.values[j] = s.values[j], s.values[i] +} + +func (s *BvhTestSotable) Less(i, j int) bool { + return s.values[i] < s.values[j] +} + +func (s *BvhTestSotable) Check(i int) bool { + return s.values[i] % 2 == 0 +} + +func TestNthElement(t *testing.T) { + numTrials := 100 + for trial := 0; trial < numTrials; trial++ { + numElems := 100 + values := make([]int, numElems) + for i := 0; i < numElems; i++ { + values[i] = rand.Intn(numElems) + } + sortable := NewBvhTestSorable(values) + splitPos := 40 //rand.Intn(numElems) + SliceNthElement(sortable, 0, numElems, splitPos) + + success := true + for i := 0; i < splitPos; i++ { + for j := splitPos; j < numElems; j++ { + if values[i] > values[j] { + t.Errorf("Condition violated: %d < %d", values[i], values[j]) + success = false + } + } + } + + if !success { + t.Errorf("%v", values) + } } +} - k := rand.Intn(3) - a := &AxisSorter{items, k} - sort.Sort(a) +func TestPartition(t *testing.T) { + numTrials := 100 + for trial := 0; trial < numTrials; trial++ { + numElems := 100 + values := make([]int, numElems) + for i := 0; i < numElems; i++ { + values[i] = rand.Intn(numElems) + } + sortable := NewBvhTestSorable(values) - for i := 0; i < len(items) - 1; i++ { - if items[i].v.NthElement(k) >= items[i + 1].v.NthElement(k) { - t.Errorf("Item sorting failed!") - break + splitPos := SlicePartition(sortable) + for i := 0; i < splitPos; i++ { + if !sortable.Check(i) { + t.Errorf("Condition violated: %d %% 2 == 0", values[i]) + } + } + for i := splitPos; i < numElems; i++ { + if sortable.Check(i) { + t.Errorf("Condition violated: %d %% 2 != 0", values[i]) + } } } } @@ -36,7 +90,7 @@ func TestBvhIntersection(t *testing.T) { triMesh := NewTriMeshFromFile("../../data/gopher.obj") bvh := NewBvh(triMesh.Primitives) - numTrials := 100 + numTrials := 1000 for trial := 0; trial < numTrials; trial++ { org := NewVector3d( rand.Float64(), diff --git a/src/accelerator/bvh_utils.go b/src/accelerator/bvh_utils.go new file mode 100644 index 0000000..59e65dd --- /dev/null +++ b/src/accelerator/bvh_utils.go @@ -0,0 +1,162 @@ +package accelerator + +import ( + "math/rand" + "sort" + . "core" +) + +const ( + BVH_SPLIT_METHOD_EQUAL_COUNT = iota + BVH_SPLIT_METHOD_SAH +) + +// ----------------------------------------------------------------------------- +// BVH node +// ----------------------------------------------------------------------------- +type BvhNode struct { + left, right *BvhNode + primId int + bbox *Bounds3d +} + +func NewLeafNode(primId int, b *Bounds3d) *BvhNode { + node := &BvhNode{} + node.left = nil + node.right = nil + node.primId = primId + node.bbox = b + return node +} + +func NewForkNode(left *BvhNode, right *BvhNode, b *Bounds3d) *BvhNode { + node := &BvhNode{} + node.left = left + node.right = right + node.primId = -1 + node.bbox = b + return node +} + +func (node *BvhNode) IsLeaf() bool { + return node.left == nil && node.right == nil +} + +// ----------------------------------------------------------------------------- +// BVH primitive info +// ----------------------------------------------------------------------------- +type BvhPrimitiveInfo struct { + primitiveId int + centroid *Vector3d + bounds *Bounds3d +} + +func NewBvhPrimitiveInfo(primitiveId int, b *Bounds3d) *BvhPrimitiveInfo { + info := &BvhPrimitiveInfo{} + info.primitiveId = primitiveId + info.bounds = b + info.centroid = b.Centroid() + return info +} + +type Partitionable interface { + sort.Interface + Check(i int) bool +} + +type BvhPrimitiveSortable struct { + items []*BvhPrimitiveInfo + axis int + checker func(*BvhPrimitiveInfo) bool +} + +func NewBvhPrimitiveSortable(items []*BvhPrimitiveInfo, axis int, checker func(*BvhPrimitiveInfo) bool) *BvhPrimitiveSortable { + s := &BvhPrimitiveSortable{} + s.items = items + s.axis = axis + s.checker = checker + return s +} + +func (s *BvhPrimitiveSortable) Len() int { + return len(s.items) +} + +func (s *BvhPrimitiveSortable) Swap(i, j int) { + s.items[i], s.items[j] = s.items[j], s.items[i] +} + +func (s *BvhPrimitiveSortable) Less(i, j int) bool { + return s.items[i].centroid.NthElement(s.axis) < s.items[j].centroid.NthElement(s.axis) +} + +func (s *BvhPrimitiveSortable) Check(i int) bool { + return s.checker(s.items[i]) +} + +// ----------------------------------------------------------------------------- +// Bucket info +// ----------------------------------------------------------------------------- +type BucketInfo struct { + count int + bounds *Bounds3d +} + +func NewBucketInfo(count int, bounds *Bounds3d) *BucketInfo { + b := &BucketInfo{} + b.count = count + b.bounds = bounds + return b +} + +// ----------------------------------------------------------------------------- +// Utility methods +// ----------------------------------------------------------------------------- +func SlicePartition(items Partitionable) int { + k := 0 + for i := 0; i < items.Len(); i++ { + if items.Check(i) { + if i != k { + items.Swap(i, k) + } + k += 1 + } + } + return k +} + +func SliceNthElement(items sort.Interface, start, end, N int) { + if end - start <= 1 || end - start <= N || N == 0 { + return + } else if end - start == 2 { + if !items.Less(start, start + 1) { + items.Swap(start, start + 1) + } + return + } + + pivot := start + rand.Intn(end - start) + if pivot != end - 1 { + items.Swap(pivot, end - 1) + pivot = end - 1 + } + + k := start + for i := start; i < end - 1; i++ { + if items.Less(i, pivot) { + if i != k { + items.Swap(i, k) + } + k += 1 + } + } + if k != pivot { + items.Swap(k, pivot) + } + + if (k - start) < N { + SliceNthElement(items, k + 1, end, N - (k - start + 1)) + } else { + SliceNthElement(items, start, k, N) + } +} diff --git a/src/core/bounds3d.go b/src/core/bounds3d.go index 4a01535..9158925 100644 --- a/src/core/bounds3d.go +++ b/src/core/bounds3d.go @@ -22,12 +22,26 @@ func NewBounds3dMinMax(minPos, maxPos *Vector3d) *Bounds3d { return b } -func (b1 *Bounds3d) Merge(b2 *Bounds3d) { - b1.MinPos = b1.MinPos.Minimum(b2.MinPos) - b1.MaxPos = b1.MaxPos.Maximum(b2.MaxPos) +func (b1 *Bounds3d) Merge(b2 *Bounds3d) *Bounds3d { + ret := &Bounds3d{} + ret.MinPos = b1.MinPos.Minimum(b2.MinPos) + ret.MaxPos = b1.MaxPos.Maximum(b2.MaxPos) + return ret } -func (b *Bounds3d) Center() *Vector3d { +func (b *Bounds3d) MergePoint(v *Vector3d) *Bounds3d { + ret := &Bounds3d{} + ret.MinPos = b.MinPos.Minimum(v) + ret.MaxPos = b.MaxPos.Maximum(v) + return ret +} + +func (b *Bounds3d) Area() Float { + diff := b.MaxPos.Subtract(b.MinPos) + return math.Abs(diff.X * diff.Y * diff.Z) +} + +func (b *Bounds3d) Centroid() *Vector3d { return b.MinPos.Add(b.MaxPos).Scale(0.5) } diff --git a/src/core/bounds3d_test.go b/src/core/bounds3d_test.go index 4247b23..a76fd2e 100644 --- a/src/core/bounds3d_test.go +++ b/src/core/bounds3d_test.go @@ -1,6 +1,7 @@ package core import ( + "math/rand" "testing" ) @@ -30,6 +31,23 @@ func TestMerge(t *testing.T) { } } +func TestMergePoint(t *testing.T) { + b := NewBounds3d() + pMin := NewVector3d(Infinity, Infinity, Infinity) + pMax := NewVector3d(-Infinity, -Infinity, -Infinity) + for i := 0; i < 100; i++ { + p := NewVector3d(rand.Float64(), rand.Float64(),rand.Float64()) + b = b.MergePoint(p) + pMin = pMin.Minimum(p) + pMax = pMax.Maximum(p) + } + + if !pMin.Equals(b.MinPos) || !pMax.Equals(b.MaxPos) { + t.Errorf("MergePoint failed: (%v, %v) != (%v, %v)", pMin, pMax, b.MinPos, b.MaxPos) + } + +} + func TestBounds3dIntersect(t *testing.T) { b1 := NewBounds3dMinMax( NewVector3d(0.0, 0.0, 0.0), @@ -50,6 +68,17 @@ func TestBounds3dIntersect(t *testing.T) { } } +func TestBounds3dCentroid(t *testing.T) { + b := NewBounds3dMinMax( + NewVector3d(1.0, 2.0, 3.0), + NewVector3d(2.0, 3.0, 4.0), + ) + c := b.Centroid() + if c.X != 1.5 || c.Y != 2.5 || c.Z != 3.5 { + t.Errorf("Unexpected centroid: %v", c) + } +} + func TestMaxExtent(t *testing.T) { b1 := NewBounds3dMinMax( NewVector3d(1.0, 2.0, 3.0), diff --git a/src/main.go b/src/main.go index ca42edf..1f6ff22 100644 --- a/src/main.go +++ b/src/main.go @@ -14,10 +14,8 @@ import ( . "sensor" . "sampler" . "integrator" - // "os" - // "log" - // "os/signal" - // "runtime/pprof" + "os/signal" + "runtime/pprof" ) var objFiles = []string{ @@ -40,22 +38,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