package assets

import "core:bytes"
import "core:encoding/csv"
import "core:io"
import "core:log"
import lg "core:math/linalg"
import "core:strconv"
import "game:debug"
import "game:halfedge"
import "game:physics/collision"
import rl "libs:raylib"
import rlgl "libs:raylib/rlgl"

CSV_Parse_Error :: enum {
	Ok,
	TooManyColumns,
	ExpectedNumber,
}

parse_csv_1d :: proc(
	data: []byte,
	allocator := context.allocator,
) -> (
	values: []f32,
	error: CSV_Parse_Error,
) {
	bytes_reader: bytes.Reader
	bytes.reader_init(&bytes_reader, data)
	bytes_stream := bytes.reader_to_stream(&bytes_reader)

	csv_reader: csv.Reader
	csv.reader_init(&csv_reader, bytes_stream, context.temp_allocator)
	defer csv.reader_destroy(&csv_reader)

	tmp_result := make([dynamic]f32, context.temp_allocator)

	skipped_header := false
	for {
		row, err := csv.read(&csv_reader, context.temp_allocator)
		if err != nil {
			if err != io.Error.EOF {
				log.warnf("Failed to read curve %v", err)
			}
			break
		}

		if len(row) != 1 {
			log.warnf("expected 1 columns, got %v", len(row))
			error = .TooManyColumns
			break
		}

		val, ok := strconv.parse_f64(row[0])
		if !ok {
			if skipped_header {
				log.warnf("Expected numbers, got %s", row[1])
				error = .ExpectedNumber
				break
			}
			skipped_header = true
			continue
		}

		append(&tmp_result, f32(val))
	}

	values = make([]f32, len(tmp_result), allocator)
	copy(values, tmp_result[:])

	return
}

parse_csv_2d :: proc(
	data: []byte,
	allocator := context.allocator,
) -> (
	curve: Curve_2D,
	error: CSV_Parse_Error,
) {
	bytes_reader: bytes.Reader
	bytes.reader_init(&bytes_reader, data)
	bytes_stream := bytes.reader_to_stream(&bytes_reader)

	csv_reader: csv.Reader
	csv.reader_init(&csv_reader, bytes_stream, context.temp_allocator)
	defer csv.reader_destroy(&csv_reader)

	tmp_result := make([dynamic][2]f32, context.temp_allocator)

	skipped_header := false
	for {
		row, err := csv.read(&csv_reader, context.temp_allocator)
		if err != nil {
			if err != io.Error.EOF {
				log.warnf("Failed to read curve %v", err)
			}
			break
		}

		if len(row) != 2 {
			log.warnf("Curve expected 2 columns, got %v", len(row))
			error = .TooManyColumns
			break
		}

		ok: bool
		key: f64
		val: f64
		key, ok = strconv.parse_f64(row[0])
		if !ok {
			if skipped_header {
				log.warnf("Curve expected numbers, got %s", row[0])
				error = .ExpectedNumber
				break
			}
			skipped_header = true
			continue
		}
		val, ok = strconv.parse_f64(row[1])
		if !ok {
			if skipped_header {
				log.warnf("Curve expected numbers, got %s", row[1])
				error = .ExpectedNumber
				break
			}
			skipped_header = true
			continue
		}

		append(&tmp_result, [2]f32{f32(key), f32(val)})
	}

	curve = make([][2]f32, len(tmp_result), allocator)
	copy(curve, tmp_result[:])

	return
}

parse_convex :: proc(bytes: []byte, allocator := context.allocator) -> (Loaded_Convex, bool) {
	Parse_Ctx :: struct {
		bytes: []byte,
		it:    int,
		line:  int,
	}

	advance :: proc(ctx: ^Parse_Ctx, by: int = 1) -> bool {
		ctx.it = min(ctx.it + by, len(ctx.bytes) + 1)
		return ctx.it < len(ctx.bytes)
	}

	is_whitespace :: proc(b: byte) -> bool {
		return b == ' ' || b == '\t' || b == '\r' || b == '\n'
	}

	skip_line :: proc(ctx: ^Parse_Ctx) {
		for ctx.it < len(ctx.bytes) && ctx.bytes[ctx.it] != '\n' {
			advance(ctx) or_break
		}
		advance(ctx)
		ctx.line += 1
	}

	skip_whitespase :: proc(ctx: ^Parse_Ctx) {
		switch ctx.bytes[ctx.it] {
		case ' ', '\t', '\r', '\n':
			if ctx.bytes[ctx.it] == '\n' {
				ctx.line += 1
			}
			advance(ctx) or_break
		case '#':
			skip_line(ctx)
		}
	}

	Edge :: [2]u16
	edges_map := make_map(map[Edge]halfedge.Edge_Index, context.temp_allocator)

	edges := make_dynamic_array([dynamic]halfedge.Half_Edge, context.temp_allocator)
	vertices := make_dynamic_array([dynamic]halfedge.Vertex, context.temp_allocator)
	faces := make_dynamic_array([dynamic]halfedge.Face, context.temp_allocator)
	min_pos, max_pos: rl.Vector3 = max(f32), min(f32)

	// Parse obj file directly into halfedge data structure
	{
		ctx := Parse_Ctx {
			bytes = bytes,
			line  = 1,
		}

		for ctx.it < len(ctx.bytes) {
			skip_whitespase(&ctx)
			switch ctx.bytes[ctx.it] {
			case 'v':
				// vertex
				advance(&ctx) or_break

				vertex: rl.Vector3

				coord_idx := 0
				for ctx.bytes[ctx.it] != '\n' && ctx.bytes[ctx.it] != '\r' {
					skip_whitespase(&ctx)
					s := string(ctx.bytes[ctx.it:])
					coord_val, nr, ok := strconv.parse_f32_prefix(s)
					if !ok {
						log.errorf(
							"failed to parse float %v %s at line %d",
							coord_idx,
							ctx.bytes[ctx.it:][:12],
							ctx.line,
						)
						return {}, false
					}
					advance(&ctx, nr) or_break

					vertex[coord_idx] = coord_val
					coord_idx += 1
				}
				append(&vertices, halfedge.Vertex{pos = vertex, edge = -1})
				min_pos = lg.min(vertex, min_pos)
				max_pos = lg.max(vertex, max_pos)

				if ctx.bytes[ctx.it] == '\r' {
					advance(&ctx)
				}
				advance(&ctx)
				ctx.line += 1
			case 'f':
				advance(&ctx) or_break

				MAX_FACE_VERTS :: 10

				indices_buf: [MAX_FACE_VERTS]u16
				index_count := 0

				for ctx.bytes[ctx.it] != '\n' && ctx.bytes[ctx.it] != '\r' {
					skip_whitespase(&ctx)
					index_f, nr, ok := strconv.parse_f32_prefix(string(ctx.bytes[ctx.it:]))
					if !ok {
						log.errorf("failed to parse index at line %d", ctx.line)
						return {}, false
					}
					advance(&ctx, nr) or_break
					index := u16(index_f) - 1
					indices_buf[index_count] = u16(index)
					index_count += 1
				}
				if ctx.bytes[ctx.it] == '\r' {
					advance(&ctx)
				}
				advance(&ctx)
				ctx.line += 1

				assert(index_count >= 3)
				indices := indices_buf[:index_count]

				append(&faces, halfedge.Face{})
				face_idx := len(faces) - 1
				face := &faces[face_idx]

				first_edge_idx := len(edges)

				face.edge = halfedge.Edge_Index(first_edge_idx)

				plane: collision.Plane
				{
					i1, i2, i3 := indices[0], indices[1], indices[2]
					v1, v2, v3 := vertices[i1].pos, vertices[i2].pos, vertices[i3].pos

					plane = collision.plane_from_point_normal(
						v1,
						lg.normalize0(lg.cross(v2 - v1, v3 - v1)),
					)
				}
				face.normal = plane.normal

				for index in indices[3:] {
					assert(
						abs(collision.signed_distance_plane(vertices[index].pos, plane)) < 0.01,
						"mesh has non planar faces",
					)
				}

				first_vert_pos := vertices[indices[0]].pos

				for i in 0 ..< len(indices) {
					edge_idx := halfedge.Edge_Index(first_edge_idx + i)
					prev_edge_relative := i == 0 ? len(indices) - 1 : i - 1
					next_edge_relative := (i + 1) % len(indices)
					i1, i2 := indices[i], indices[next_edge_relative]
					v1, v2 := &vertices[i1], &vertices[i2]

					assert(
						lg.dot(
							lg.cross(v1.pos - first_vert_pos, v2.pos - first_vert_pos),
							plane.normal,
						) >=
						0,
						"non convex face or non ccw winding",
					)

					if v1.edge == -1 {
						v1.edge = edge_idx
					}

					edge := halfedge.Half_Edge {
						origin = halfedge.Vertex_Index(i1),
						face   = halfedge.Face_Index(face_idx),
						twin   = -1,
						next   = halfedge.Edge_Index(first_edge_idx + next_edge_relative),
						prev   = halfedge.Edge_Index(first_edge_idx + prev_edge_relative),
					}

					stable_index := [2]u16{min(i1, i2), max(i1, i2)}
					if stable_index in edges_map {
						edge.twin = edges_map[stable_index]
						twin_edge := &edges[edge.twin]
						assert(twin_edge.twin == -1, "edge has more than two faces attached")
						twin_edge.twin = edge_idx
					} else {
						edges_map[stable_index] = edge_idx
					}

					append(&edges, edge)
				}
			case:
				skip_line(&ctx)
			}
		}
	}

	center := (max_pos + min_pos) * 0.5
	extent := (max_pos - min_pos) * 0.5

	center_of_mass: rl.Vector3

	final_vertices := make([]halfedge.Vertex, len(vertices), allocator)
	final_edges := make([]halfedge.Half_Edge, len(edges), allocator)
	final_faces := make([]halfedge.Face, len(faces), allocator)
	copy(final_vertices, vertices[:])
	copy(final_edges, edges[:])
	copy(final_faces, faces[:])

	mesh := halfedge.Half_Edge_Mesh {
		vertices = final_vertices,
		edges    = final_edges,
		faces    = final_faces,
		center   = center,
		extent   = extent,
	}

	// Center of mass calculation
	total_volume := f32(0.0)
	{
		tri_idx := 0
		for face_idx in 0 ..< len(faces) {
			face := faces[face_idx]
			// for all triangles
			it := halfedge.iterator_face_edges(mesh, halfedge.Face_Index(face_idx))
			i := 0
			tri: [3]rl.Vector3
			for edge in halfedge.iterate_next_edge(&it) {
				switch i {
				case 0 ..< 3:
					tri[i] = mesh.vertices[edge.origin].pos
				case:
					tri[1] = tri[2]
					tri[2] = mesh.vertices[edge.origin].pos
				}

				if i >= 2 {
					plane := collision.plane_from_point_normal(tri[0], -face.normal)

					h := max(0, collision.signed_distance_plane(center, plane))
					tri_area :=
						lg.dot(lg.cross(tri[1] - tri[0], tri[2] - tri[0]), face.normal) * 0.5
					tetra_volume := 1.0 / 3.0 * tri_area * h
					total_volume += tetra_volume

					tetra_centroid := (tri[0] + tri[1] + tri[2] + center) * 0.25
					center_of_mass += tetra_volume * tetra_centroid

					tri_idx += 1
				}

				i += 1
			}
		}
	}

	assert(total_volume > 0, "degenerate convex hull")
	center_of_mass /= total_volume

	inertia_tensor: lg.Matrix3f32
	// Find inertia tensor
	{
		tri_idx := 0
		for face_idx in 0 ..< len(faces) {
			// for all triangles
			it := halfedge.iterator_face_edges(mesh, halfedge.Face_Index(face_idx))
			i := 0
			tri: [3]rl.Vector3
			for edge in halfedge.iterate_next_edge(&it) {
				switch i {
				case 0 ..< 3:
					tri[i] = mesh.vertices[edge.origin].pos
				case:
					tri[1] = tri[2]
					tri[2] = mesh.vertices[edge.origin].pos
				}

				if i >= 2 {
					tet := Tetrahedron {
						p = {tri[0], tri[1], tri[2], center_of_mass},
					}

					inertia_tensor += tetrahedron_inertia_tensor(tet, center_of_mass)

					tri_idx += 1
				}

				i += 1
			}
		}
	}
	inertia_tensor = inertia_tensor

	return Loaded_Convex {
			mesh = mesh,
			center_of_mass = center_of_mass,
			inertia_tensor = inertia_tensor,
			total_volume = total_volume,
		},
		true
}

// TODO: move convex stuff out of assets.odin
Tetrahedron :: struct {
	p: [4]rl.Vector3,
}

tetrahedron_volume :: #force_inline proc(tet: Tetrahedron) -> f32 {
	return(
		1.0 /
		6.0 *
		abs(lg.dot(lg.cross(tet.p[1] - tet.p[0], tet.p[2] - tet.p[0]), tet.p[3] - tet.p[0])) \
	)
}

square :: #force_inline proc(val: f32) -> f32 {
	return val * val
}

tetrahedron_inertia_tensor :: proc(tet: Tetrahedron, o: rl.Vector3) -> lg.Matrix3f32 {
	p1, p2, p3, p4 := tet.p[0] - o, tet.p[1] - o, tet.p[2] - o, tet.p[3] - o
	// Jacobian determinant is 6*Volume
	det_j := abs(6.0 * tetrahedron_volume(tet))

	moment_of_inertia_term :: proc(p1, p2, p3, p4: rl.Vector3, axis: int) -> f32 {
		return(
			square(p1[axis]) +
			p1[axis] * p2[axis] +
			square(p2[axis]) +
			p1[axis] * p3[axis] +
			p2[axis] * p3[axis] +
			square(p3[axis]) +
			p1[axis] * p4[axis] +
			p2[axis] * p4[axis] +
			p3[axis] * p4[axis] +
			square(p4[axis]) \
		)
	}

	product_of_inertia_term :: proc(p1, p2, p3, p4: rl.Vector3, axis1, axis2: int) -> f32 {
		return(
			2.0 * p1[axis1] * p1[axis2] +
			p2[axis1] * p1[axis2] +
			p3[axis1] * p1[axis2] +
			p4[axis1] * p1[axis2] +
			p1[axis1] * p2[axis2] +
			2.0 * p2[axis1] * p2[axis2] +
			p3[axis1] * p2[axis2] +
			p4[axis1] * p2[axis2] +
			p1[axis1] * p3[axis2] +
			p2[axis1] * p3[axis2] +
			2.0 * p3[axis1] * p3[axis2] +
			p4[axis1] * p3[axis2] +
			p1[axis1] * p4[axis2] +
			p2[axis1] * p4[axis2] +
			p3[axis1] * p4[axis2] +
			2.0 * p4[axis1] * p4[axis2] \
		)
	}

	MOMENT_OF_INERTIA_DENOM :: 1.0 / 60.0
	PRODUCT_OF_INERTIA_DENOM :: 1.0 / 120.0

	x_term := moment_of_inertia_term(p1, p2, p3, p4, 0)
	y_term := moment_of_inertia_term(p1, p2, p3, p4, 1)
	z_term := moment_of_inertia_term(p1, p2, p3, p4, 2)

	// Moments of intertia with respect to XYZ
	// Integral(y^2 + z^2)
	a := det_j * (y_term + z_term) * MOMENT_OF_INERTIA_DENOM
	// Integral(x^2 + z^2)
	b := det_j * (x_term + z_term) * MOMENT_OF_INERTIA_DENOM
	// Integral(x^2 + y^2)
	c := det_j * (x_term + y_term) * MOMENT_OF_INERTIA_DENOM

	// Products of inertia
	a_ := product_of_inertia_term(p1, p2, p3, p4, axis1 = 1, axis2 = 2) * PRODUCT_OF_INERTIA_DENOM
	b_ := product_of_inertia_term(p1, p2, p3, p4, axis1 = 0, axis2 = 2) * PRODUCT_OF_INERTIA_DENOM
	c_ := product_of_inertia_term(p1, p2, p3, p4, axis1 = 0, axis2 = 1) * PRODUCT_OF_INERTIA_DENOM

	return {a, -b_, -c_, -b_, b, -a_, -c_, -a_, c}
}

debug_draw_tetrahedron_wires :: proc(tri: [3]rl.Vector3, p: rl.Vector3, color: rl.Color) {
	rlgl.Begin(rlgl.LINES)
	defer rlgl.End()

	debug.rlgl_color(color)

	debug.rlgl_vertex3v2(tri[0], tri[1])
	debug.rlgl_vertex3v2(tri[1], tri[2])
	debug.rlgl_vertex3v2(tri[2], tri[0])
	debug.rlgl_vertex3v2(tri[0], p)
	debug.rlgl_vertex3v2(tri[1], p)
	debug.rlgl_vertex3v2(tri[2], p)
}