view cmd/octree2contour/vertex.go @ 687:be90ab542aa7 octree

octree: contouring: Do the math to calculate the intersection points of the triangles and the planes.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 20 Sep 2018 11:32:03 +0200
parents 120a82bd9953
children 614135d69823
line wrap: on
line source

package main

type vertex struct {
	x float64
	y float64
	z float64
}

func (v *vertex) minimize(w vertex) {
	if w.x < v.x {
		v.x = w.x
	}
	if w.y < v.y {
		v.y = w.y
	}
	if w.z < v.z {
		v.z = w.z
	}
}

func (v *vertex) maximize(w vertex) {
	if w.x > v.x {
		v.x = w.x
	}
	if w.y > v.y {
		v.y = w.y
	}
	if w.z > v.z {
		v.z = w.z
	}
}

func (v vertex) sub(w vertex) vertex {
	return vertex{
		v.x - w.x,
		v.y - w.y,
		v.z - w.z,
	}
}

func (v vertex) add(w vertex) vertex {
	return vertex{
		v.x + w.x,
		v.y + w.y,
		v.z + w.z,
	}
}

func (v vertex) scale(s float64) vertex {
	return vertex{
		s * v.x,
		s * v.y,
		s * v.z,
	}
}

func interpolate(v1, v2 vertex) func(vertex) vertex {
	v2 = v2.sub(v1)
	return func(s vertex) vertex {
		return vertex{
			v2.x*s.x + v1.x,
			v2.y*s.y + v1.y,
			v2.z*s.z + v1.z,
		}
	}
}

func (a vertex) less(b vertex) bool {
	return a.x < b.x || a.y < b.y || a.z < b.z
}

type line [2]vertex

func newLine(p1, p2 vertex) line {
	return line{
		p2.sub(p1),
		p1,
	}
}

func (l line) eval(t float64) vertex {
	return l[0].scale(t).add(l[1])
}

func (l line) intersectH(h float64) vertex {
	t := (h - l[1].z) / l[0].z
	return l.eval(t)
}

type triangle [3]vertex

func side(z, h float64) int {
	switch {
	case z < h:
		return -1
	case z > h:
		return +1
	}
	return 0
}

func (t *triangle) intersectH(h float64, points []vertex) []vertex {
	sides := [3]int{
		side(t[0].z, h),
		side(t[1].z, h),
		side(t[2].z, h),
	}

	for i := 0; i < 3; i++ {
		j := (i + 1) % 3
		si := sides[i]
		sj := sides[j]

		switch {
		case si == 0 && sj == 0:
			// both on plane
			points = append(points, t[i], t[j])
		case si == 0 && sj != 0:
			// first on plane
			points = append(points, t[i])
		case si != 0 && sj == 0:
			// second on plane
			points = append(points, t[j])
		case si == sj:
			// both on same side
		default:
			// real intersection
			v := newLine(t[i], t[j]).intersectH(h)
			points = append(points, v)
		}
	}

	return points
}