simple rasterizer

raster.go

package main

import (
	"fmt"
	"image"
	"image/color"
	"image/png"
	"log"
	"math"
	"net/http"
)

type Vec3 struct {
	X float64
	Y float64
	Z float64
}

type Vec4 struct {
	X float64
	Y float64
	Z float64
	W float64
}

type Mat4 [16]float64

type Attributes struct {
	Position Vec3
	Normal   Vec3
}

type Interpolated [3]float64

var cubeAttributes = []Attributes{
	// front
	{Vec3{-0.5, -0.5, -0.5}, Vec3{0, 0, -1}},
	{Vec3{0.5, -0.5, -0.5}, Vec3{0, 0, -1}},
	{Vec3{-0.5, 0.5, -0.5}, Vec3{0, 0, -1}},

	{Vec3{0.5, 0.5, -0.5}, Vec3{0, 0, -1}},
	{Vec3{0.5, -0.5, -0.5}, Vec3{0, 0, -1}},
	{Vec3{-0.5, 0.5, -0.5}, Vec3{0, 0, -1}},

	// back
	{Vec3{-0.5, -0.5, 0.5}, Vec3{0, 0, 1}},
	{Vec3{0.5, -0.5, 0.5}, Vec3{0, 0, 1}},
	{Vec3{-0.5, 0.5, 0.5}, Vec3{0, 0, 1}},

	{Vec3{0.5, 0.5, 0.5}, Vec3{0, 0, 1}},
	{Vec3{0.5, -0.5, 0.5}, Vec3{0, 0, 1}},
	{Vec3{-0.5, 0.5, 0.5}, Vec3{0, 0, 1}},

	// top
	{Vec3{-0.5, -0.5, -0.5}, Vec3{0, -1, 0}},
	{Vec3{0.5, -0.5, -0.5}, Vec3{0, -1, 0}},
	{Vec3{-0.5, -0.5, 0.5}, Vec3{0, -1, 0}},

	{Vec3{0.5, -0.5, 0.5}, Vec3{0, -1, 0}},
	{Vec3{0.5, -0.5, -0.5}, Vec3{0, -1, 0}},
	{Vec3{-0.5, -0.5, 0.5}, Vec3{0, -1, 0}},

	// bottom
	{Vec3{-0.5, 0.5, -0.5}, Vec3{0, 1, 0}},
	{Vec3{0.5, 0.5, -0.5}, Vec3{0, 1, 0}},
	{Vec3{-0.5, 0.5, 0.5}, Vec3{0, 1, 0}},

	{Vec3{0.5, 0.5, 0.5}, Vec3{0, 1, 0}},
	{Vec3{0.5, 0.5, -0.5}, Vec3{0, 1, 0}},
	{Vec3{-0.5, 0.5, 0.5}, Vec3{0, 1, 0}},

	// left
	{Vec3{-0.5, -0.5, -0.5}, Vec3{-1, 0, 0}},
	{Vec3{-0.5, -0.5, 0.5}, Vec3{-1, 0, 0}},
	{Vec3{-0.5, 0.5, -0.5}, Vec3{-1, 0, 0}},

	{Vec3{-0.5, 0.5, 0.5}, Vec3{-1, 0, 0}},
	{Vec3{-0.5, -0.5, 0.5}, Vec3{-1, 0, 0}},
	{Vec3{-0.5, 0.5, -0.5}, Vec3{-1, 0, 0}},

	// right
	{Vec3{0.5, -0.5, -0.5}, Vec3{1, 0, 0}},
	{Vec3{0.5, -0.5, 0.5}, Vec3{1, 0, 0}},
	{Vec3{0.5, 0.5, -0.5}, Vec3{1, 0, 0}},

	{Vec3{0.5, 0.5, 0.5}, Vec3{1, 0, 0}},
	{Vec3{0.5, -0.5, 0.5}, Vec3{1, 0, 0}},
	{Vec3{0.5, 0.5, -0.5}, Vec3{1, 0, 0}},
}

func (a Vec3) Dot(b Vec3) float64 {
	return a.X*b.X + a.Y*b.Y + a.Z*b.Z
}

func (v Vec3) MulScalar(f float64) Vec3 {
	return Vec3{v.X * f, v.Y * f, v.Z * f}
}

func (v Vec3) Normalize() Vec3 {
	mag := math.Sqrt(v.X*v.X + v.Y*v.Y + v.Z*v.Z)
	return v.MulScalar(1 / mag)
}

func (v Vec4) Homogenize() Vec4 {
	return Vec4{v.X / v.W, v.Y / v.W, v.Z / v.W, 1}
}

func (a Mat4) MulMat4(b Mat4) Mat4 {
	return Mat4{
		a[0]*b[0] + a[1]*b[4] + a[2]*b[8] + a[3]*b[12],
		a[0]*b[1] + a[1]*b[5] + a[2]*b[9] + a[3]*b[13],
		a[0]*b[2] + a[1]*b[6] + a[2]*b[10] + a[3]*b[14],
		a[0]*b[3] + a[1]*b[7] + a[2]*b[11] + a[3]*b[15],
		a[4]*b[0] + a[5]*b[4] + a[6]*b[8] + a[7]*b[12],
		a[4]*b[1] + a[5]*b[5] + a[6]*b[9] + a[7]*b[13],
		a[4]*b[2] + a[5]*b[6] + a[6]*b[10] + a[7]*b[14],
		a[4]*b[3] + a[5]*b[7] + a[6]*b[11] + a[7]*b[15],
		a[8]*b[0] + a[9]*b[4] + a[10]*b[8] + a[11]*b[12],
		a[8]*b[1] + a[9]*b[5] + a[10]*b[9] + a[11]*b[13],
		a[8]*b[2] + a[9]*b[6] + a[10]*b[10] + a[11]*b[14],
		a[8]*b[3] + a[9]*b[7] + a[10]*b[11] + a[11]*b[15],
		a[12]*b[0] + a[13]*b[4] + a[14]*b[8] + a[15]*b[12],
		a[12]*b[1] + a[13]*b[5] + a[14]*b[9] + a[15]*b[13],
		a[12]*b[2] + a[13]*b[6] + a[14]*b[10] + a[15]*b[14],
		a[12]*b[3] + a[13]*b[7] + a[14]*b[11] + a[15]*b[15],
	}
}

func (m Mat4) MulVec4(v Vec4) Vec4 {
	return Vec4{
		v.X*m[0] + v.Y*m[1] + v.Z*m[2] + v.W*m[3],
		v.X*m[4] + v.Y*m[5] + v.Z*m[6] + v.W*m[7],
		v.X*m[8] + v.Y*m[9] + v.Z*m[10] + v.W*m[11],
		v.X*m[12] + v.Y*m[13] + v.Z*m[14] + v.W*m[15],
	}
}

func MakeYRotation(radians float64) Mat4 {
	return Mat4{
		math.Cos(radians), 0, math.Sin(radians), 0,
		0, 1, 0, 0,
		-math.Sin(radians), 0, math.Cos(radians), 0,
		0, 0, 0, 1,
	}
}

func MakeTranslation(x, y, z float64) Mat4 {
	return Mat4{
		1, 0, 0, x,
		0, 1, 0, y,
		0, 0, 1, z,
		0, 0, 0, 1,
	}
}

func MakeScale(x, y, z float64) Mat4 {
	return Mat4{
		x, 0, 0, 0,
		0, y, 0, 0,
		0, 0, z, 0,
		0, 0, 0, 1,
	}
}

func MakeOrtho(left, right, bottom, top, near, far float64) Mat4 {
	return MakeScale(2/(right-left), 2/(top-bottom), 2/(far-near)).MulMat4(MakeTranslation(-(right+left)/2, -(top+bottom)/2, -(far+near)/2))
}

func MakePerspective(fovDegrees, near, far float64) Mat4 {
	fovRadians := fovDegrees * math.Pi / 180
	scale := 1 / math.Tan(fovRadians/2)
	return Mat4{
		scale, 0, 0, 0,
		0, scale, 0, 0,
		0, 0, (far + near) / (far - near), 2 * near * far / (near - far),
		0, 0, 1, 0,
	}
}

func render(img *image.NRGBA, vertexShader func(a Attributes) (Vec4, Interpolated), fragmentShader func(pos Vec4, interp Interpolated) color.Color) error {
	size := img.Bounds().Max.X

	// clear the image
	for py := 0; py < size; py++ {
		for px := 0; px < size; px++ {
			img.Set(px, py, color.NRGBA{127, 127, 127, 255})
		}
	}

	vertices := []Vec4{}
	interps := []Interpolated{}
	for _, v := range cubeAttributes {
		vertex, interp := vertexShader(v)
		vertices = append(vertices, vertex.Homogenize())
		interps = append(interps, interp)
	}

	// depth buffer so that we can draw triangles in any order and they don't overlap incorrectly
	depth := make([]float64, size*size)
	for i := range depth {
		depth[i] = 1
	}

	// find which side of a line a point is on using magnitude of cross product
	side := func(x0, y0, x1, y1, px, py float64) bool {
		return (x1-x0)*(py-y0)-(y1-y0)*(px-x0) > 0
	}

	// generate fragments
	for i := 0; i < len(vertices); i += 3 {
		// for _, i := range indices {
		a := vertices[i]
		b := vertices[i+1]
		c := vertices[i+2]

		interpA := interps[i]
		interpB := interps[i+1]
		interpC := interps[i+2]

		// create bounding boxes for triangles
		// it's really slow without bounding boxes
		minX := math.Min(a.X, math.Min(b.X, c.X))
		minY := math.Min(a.Y, math.Min(b.Y, c.Y))
		maxX := math.Max(a.X, math.Max(b.X, c.X))
		maxY := math.Max(a.Y, math.Max(b.Y, c.Y))

		minPx := int(math.Floor((minX+1)/2*float64(size) - 0.5))
		if minPx < 0 {
			minPx = 0
		}
		minPy := int(math.Floor((minY+1)/2*float64(size) - 0.5))
		if minPy < 0 {
			minPy = 0
		}
		maxPx := int(math.Ceil((maxX+1)/2*float64(size) - 0.5))
		if maxPx >= size {
			maxPx = size - 1
		}
		maxPy := int(math.Ceil((maxY+1)/2*float64(size) - 0.5))
		if maxPy >= size {
			maxPy = size - 1
		}

		// generate all pixels that fall into this box
		for py := minPy; py < maxPy; py++ {
			for px := minPx; px < maxPx; px++ {
				// check which pixels have their center inside the triangle
				x := (float64(px)+0.5)/float64(size)*2 - 1
				y := (float64(py)+0.5)/float64(size)*2 - 1

				s0 := side(a.X, a.Y, b.X, b.Y, x, y)
				s1 := side(b.X, b.Y, c.X, c.Y, x, y)
				s2 := side(c.X, c.Y, a.X, a.Y, x, y)

				// if point p is on the same side of all line segments, it's inside the triangle
				if (s0 == s1) && (s1 == s2) {
					// calculate barycentric coordinates
					// http://en.wikipedia.org/wiki/Barycentric_coordinate_system
					denom := (b.Y-c.Y)*(a.X-c.X) + (c.X-b.X)*(a.Y-c.Y)
					b0 := ((b.Y-c.Y)*(x-c.X) + (c.X-b.X)*(y-c.Y)) / denom
					b1 := ((c.Y-a.Y)*(x-c.X) + (a.X-c.X)*(y-c.Y)) / denom
					bary := Vec3{b0, b1, 1 - b0 - b1}

					// calculate depth at x,y on the surface of the triangle
					// is this correct?
					pointDepth := bary.X*a.Z + bary.Y*b.Z + bary.Z*c.Z
					if pointDepth < depth[py*size+px] {
						depth[py*size+px] = pointDepth

						// interpolate
						interp := Interpolated{}
						for i := range interp {
							interp[i] = interpA[i]*bary.X + interpB[i]*bary.Y + interpC[i]*bary.Z
						}

						c := fragmentShader(Vec4{x, y, pointDepth, 1}, interp)
						img.Set(px, py, c)
					}
				}
			}
		}
	}

	return nil
}

func streamHandler(w http.ResponseWriter, r *http.Request) {
	const boundary = "MixedReplaceBoundary"

	w.Header().Set("Content-Type", "multipart/x-mixed-replace;boundary="+boundary)

	img := image.NewNRGBA(image.Rect(0, 0, 512, 512))

	frame := 0

	// process the vertex data
	projection := MakePerspective(90, 5, 15)
	// projection := MakeOrtho(-5, 5, -5, 5, 5, 15)
	view := MakeTranslation(0, 0, 10)
	model := MakeScale(5, 5, 5).MulMat4(MakeYRotation(float64(frame) / 10))
	transform := projection.MulMat4(view).MulMat4(model)
	normalTransform := MakeYRotation(float64(frame) / 10)

	vertexShader := func(a Attributes) (Vec4, Interpolated) {
		// calculate position
		v := transform.MulVec4(Vec4{a.Position.X, a.Position.Y, a.Position.Z, 1})

		// calculate color (interpolated across triangle)
		eye4 := normalTransform.MulVec4(Vec4{a.Normal.X, a.Normal.Y, a.Normal.Z, 1})
		eye := Vec3{eye4.X, eye4.Y, eye4.Z}
		light := Vec3{-1, -1, -1}
		dotProduct := math.Max(0, eye.Normalize().Dot(light.Normalize()))
		diffuseColor := Vec3{0.4, 0.4, 1.0}
		c := diffuseColor.MulScalar(dotProduct)
		interp := Interpolated{c.X, c.Y, c.Z}
		return v, interp
	}

	fragmentShader := func(pos Vec4, interp Interpolated) color.Color {
		return color.NRGBA{
			uint8(interp[0] * 255),
			uint8(interp[1] * 255),
			uint8(interp[2] * 255),
			255,
		}
	}

	// render loop
	for {
		model = MakeScale(5, 5, 5).MulMat4(MakeYRotation(float64(frame) / 10))
		transform = projection.MulMat4(view).MulMat4(model)
		normalTransform = MakeYRotation(float64(frame) / 10)

		if err := render(img, vertexShader, fragmentShader); err != nil {
			log.Fatal(err)
		}

		w.Write([]byte(fmt.Sprintf("Content-Type: image/png\r\n\r\n")))
		encoder := png.Encoder{png.NoCompression}
		if err := encoder.Encode(w, img); err != nil {
			return
		}
		w.Write([]byte("\r\n--" + boundary + "\r\n"))
		if f, ok := w.(http.Flusher); ok {
			f.Flush()
		}

		frame++
	}
}

func indexHandler(w http.ResponseWriter, r *http.Request) {
	w.Header().Set("Content-Type", "text/html")
	w.Write([]byte(`<div style="text-align: center"><img src="/stream"/></div>`))
}

func main() {
	http.HandleFunc("/stream", streamHandler)
	http.HandleFunc("/", indexHandler)
	log.Fatal(http.ListenAndServe("127.0.0.1:8080", nil))
}

xbary.go

// show barycentric coords
c := color.NRGBA{
	uint8(bary.X * 255),
	uint8(bary.Y * 255),
	uint8(bary.Z * 255),
	255,
}

xdepth.go

// show depth buffer
for i, d := range depth {
	px := i % size
	py := i / size
	a := uint8((d + 1) / 2 * 255)
	img.Set(px, py, color.NRGBA{a, a, a, 255})
}