Heightmap to Isolines

isoline.go

package main

import (
	"fmt"
	"image"
	"image/draw"
	"image/png"
	"math"
	"os"

	"github.com/fogleman/fauxgl"
	"github.com/fogleman/gg"
)

func LoadImage(path string) (image.Image, error) {
	file, err := os.Open(path)
	if err != nil {
		return nil, err
	}
	defer file.Close()
	im, _, err := image.Decode(file)
	return im, err
}

func SavePNG(path string, im image.Image) error {
	file, err := os.Create(path)
	if err != nil {
		return err
	}
	defer file.Close()
	return png.Encode(file, im)
}

func EnsureGray(im image.Image) (*image.Gray, bool) {
	switch im := im.(type) {
	case *image.Gray:
		return im, true
	default:
		dst := image.NewGray(im.Bounds())
		draw.Draw(dst, im.Bounds(), im, image.ZP, draw.Src)
		return dst, false
	}
}

func intersectSegment(z float64, v0, v1 fauxgl.Vector) (fauxgl.Vector, bool) {
	if v0.Z == v1.Z {
		return fauxgl.Vector{}, false
	}
	t := (z - v0.Z) / (v1.Z - v0.Z)
	if t < 0 || t > 1 {
		return fauxgl.Vector{}, false
	}
	v := v0.Add(v1.Sub(v0).MulScalar(t))
	return v, true
}

func intersectTriangle(z float64, pt1, pt2, pt3 fauxgl.Vector) (fauxgl.Vector, fauxgl.Vector, bool) {
	v1, ok1 := intersectSegment(z, pt1, pt2)
	v2, ok2 := intersectSegment(z, pt2, pt3)
	v3, ok3 := intersectSegment(z, pt3, pt1)
	var p1, p2 fauxgl.Vector
	if ok1 && ok2 {
		p1, p2 = v1, v2
	} else if ok1 && ok3 {
		p1, p2 = v1, v3
	} else if ok2 && ok3 {
		p1, p2 = v2, v3
	} else {
		return fauxgl.Vector{}, fauxgl.Vector{}, false
	}
	n := fauxgl.Vector{p1.Y - p2.Y, p2.X - p1.X, 0}
	e1 := pt2.Sub(pt1)
	e2 := pt3.Sub(pt1)
	tn := e1.Cross(e2).Normalize()
	if n.Dot(tn) < 0 {
		return p1, p2, true
	} else {
		return p2, p1, true
	}
}

type Point struct {
	X, Y float64
}

type Pair struct {
	A, B Point
}

func Slice(im *image.Gray, z float64) []Pair {
	var result []Pair
	w := im.Bounds().Size().X
	h := im.Bounds().Size().Y
	for y := 0; y < h-1; y++ {
		i0 := im.PixOffset(0, y)
		i1 := im.PixOffset(0, y+1)
		for x := 1; x < w-1; x++ {
			z0 := float64(im.Pix[i0+x])
			z1 := float64(im.Pix[i0+x+1])
			z2 := float64(im.Pix[i1+x-1])
			z3 := float64(im.Pix[i1+x])
			if z0 < z && z1 < z && z2 < z && z3 < z {
				continue
			}
			if z0 > z && z1 > z && z2 > z && z3 > z {
				continue
			}
			v0 := fauxgl.Vector{float64(x), float64(y), z0}
			v1 := fauxgl.Vector{float64(x + 1), float64(y), z1}
			v2 := fauxgl.Vector{float64(x - 1), float64(y + 1), z2}
			v3 := fauxgl.Vector{float64(x), float64(y + 1), z3}
			if p1, p2, ok := intersectTriangle(z, v0, v2, v3); ok {
				pair := Pair{Point{p1.X, p1.Y}, Point{p2.X, p2.Y}}
				result = append(result, pair)
			}
			if p1, p2, ok := intersectTriangle(z, v0, v3, v1); ok {
				pair := Pair{Point{p1.X, p1.Y}, Point{p2.X, p2.Y}}
				result = append(result, pair)
			}
		}
	}
	return result
}

// func JoinPaths(paths []Path) []Path {
// 	lookup := make(map[fauxgl.Vector]Path, len(paths))
// 	for _, path := range paths {
// 		lookup[path[0]] = path
// 	}
// 	var result []Path
// 	for len(lookup) > 0 {
// 		var v fauxgl.Vector
// 		for v = range lookup {
// 			break
// 		}
// 		var path Path
// 		for {
// 			path = append(path, v)
// 			if p, ok := lookup[v]; ok {
// 				delete(lookup, v)
// 				v = p[len(p)-1]
// 			} else {
// 				break
// 			}
// 		}
// 		result = append(result, path)
// 	}
// 	return result
// }

func Render(w, h int, scale float64, pairs []Pair) image.Image {
	dc := gg.NewContext(int(float64(w)*scale), int(float64(h)*scale))
	dc.SetRGB(1, 1, 1)
	dc.Clear()
	dc.Scale(scale, scale)
	for _, pair := range pairs {
		a := pair.A
		b := pair.B
		dc.DrawLine(a.X, a.Y, b.X, b.Y)
	}
	dc.SetRGB(0, 0, 0)
	dc.Stroke()
	return dc.Image()
}

func main() {
	args := os.Args[1:]
	im, err := LoadImage(args[0])
	if err != nil {
		panic(err)
	}
	gray, _ := EnsureGray(im)
	fmt.Println(gray.Bounds())
	w := gray.Bounds().Size().X
	h := gray.Bounds().Size().Y
	var pairs []Pair
	var z0, z1 float64
	z0 = -1
	z1 = 256
	n := 16
	for i := 0; i < n; i++ {
		pct := float64(i) / float64(n-1)
		pct = math.Pow(pct, 2)
		z := z0 + (z1-z0)*pct + 1e-3
		p := Slice(gray, float64(z)+0.5)
		pairs = append(pairs, p...)
		fmt.Println(z, len(p))
	}
	iso := Render(w, h, 1, pairs)
	gg.SavePNG("out.png", iso)
}