package font

import "core:mem"
import "core:math"

// ============================================================================
// SDF (SIGNED DISTANCE FIELD) RENDERING (Phase 6)
// ============================================================================

@(private)
cuberoot :: proc(x: f32) -> f32 {
	if x < 0 {
		return -math.pow(-x, 1.0 / 3.0)
	} else {
		return math.pow(x, 1.0 / 3.0)
	}
}

// Solve cubic equation: x^3 + a*x^2 + b*x + c = 0
// Returns number of roots (1 or 3), stores roots in r
@(private)
solve_cubic :: proc(a, b, c: f32, r: [^]f32) -> i32 {
	s := -a / 3
	p := b - a * a / 3
	q := a * (2 * a * a - 9 * b) / 27 + c
	p3 := p * p * p
	d := q * q + 4 * p3 / 27
	if d >= 0 {
		z := math.sqrt(d)
		u := (-q + z) / 2
		v := (-q - z) / 2
		u = cuberoot(u)
		v = cuberoot(v)
		r[0] = s + u + v
		return 1
	} else {
		u := math.sqrt(-p / 3)
		v := math.acos(-math.sqrt(-27 / p3) * q / 2) / 3 // p3 must be negative since d is negative
		m := math.cos(v)
		n := math.cos(v - 3.141592 / 2) * 1.732050808
		r[0] = s + u * 2 * m
		r[1] = s - u * (m + n)
		r[2] = s - u * (m - n)
		return 3
	}
}

// Find intersections between a ray and a quadratic Bezier curve
// Returns number of hits (0, 1, or 2)
@(private)
ray_intersect_bezier :: proc(orig: [2]f32, ray: [2]f32, q0: [2]f32, q1: [2]f32, q2: [2]f32, hits: ^[2][2]f32) -> i32 {
	q0perp := q0[1] * ray[0] - q0[0] * ray[1]
	q1perp := q1[1] * ray[0] - q1[0] * ray[1]
	q2perp := q2[1] * ray[0] - q2[0] * ray[1]
	roperp := orig[1] * ray[0] - orig[0] * ray[1]

	a := q0perp - 2 * q1perp + q2perp
	b := q1perp - q0perp
	c := q0perp - roperp

	s0: f32 = 0
	s1: f32 = 0
	num_s: i32 = 0

	if a != 0.0 {
		discr := b * b - a * c
		if discr > 0.0 {
			rcpna := -1 / a
			d := math.sqrt(discr)
			s0 = (b + d) * rcpna
			s1 = (b - d) * rcpna
			if s0 >= 0.0 && s0 <= 1.0 {
				num_s = 1
			}
			if d > 0.0 && s1 >= 0.0 && s1 <= 1.0 {
				if num_s == 0 {
					s0 = s1
				}
				num_s += 1
			}
		}
	} else {
		// 2*b*s + c = 0
		// s = -c / (2*b)
		s0 = c / (-2 * b)
		if s0 >= 0.0 && s0 <= 1.0 {
			num_s = 1
		}
	}

	if num_s == 0 {
		return 0
	} else {
		rcp_len2 := 1 / (ray[0] * ray[0] + ray[1] * ray[1])
		rayn_x := ray[0] * rcp_len2
		rayn_y := ray[1] * rcp_len2

		q0d := q0[0] * rayn_x + q0[1] * rayn_y
		q1d := q1[0] * rayn_x + q1[1] * rayn_y
		q2d := q2[0] * rayn_x + q2[1] * rayn_y
		rod := orig[0] * rayn_x + orig[1] * rayn_y

		q10d := q1d - q0d
		q20d := q2d - q0d
		q0rd := q0d - rod

		hits[0][0] = q0rd + s0 * (2.0 - 2.0 * s0) * q10d + s0 * s0 * q20d
		hits[0][1] = a * s0 + b

		if num_s > 1 {
			hits[1][0] = q0rd + s1 * (2.0 - 2.0 * s1) * q10d + s1 * s1 * q20d
			hits[1][1] = a * s1 + b
			return 2
		} else {
			return 1
		}
	}
}

@(private)
point_equal :: proc(a: [2]f32, b: [2]f32) -> bool {
	return a[0] == b[0] && a[1] == b[1]
}

// Compute winding number by counting edge crossings from (-infinity, y) to (x, y)
@(private)
compute_crossings_x :: proc(x, y_in: f32, nverts: i32, verts: [^]Vertex) -> i32 {
	ray: [2]f32 = {1, 0}
	winding: i32 = 0

	// Make sure y never passes through a vertex of the shape
	y := y_in
	y_frac := math.mod(y, 1.0)
	if y_frac < 0.01 {
		y += 0.01
	} else if y_frac > 0.99 {
		y -= 0.01
	}

	orig: [2]f32 = {x, y}

	// Test a ray from (-infinity, y) to (x, y)
	#no_bounds_check for i: i32 = 0; i < nverts; i += 1 {
		if verts[i].type == VLINE {
			x0 := i32(verts[i - 1].x)
			y0 := i32(verts[i - 1].y)
			x1 := i32(verts[i].x)
			y1 := i32(verts[i].y)
			if y > f32(min(y0, y1)) && y < f32(max(y0, y1)) && x > f32(min(x0, x1)) {
				x_inter := (y - f32(y0)) / f32(y1 - y0) * f32(x1 - x0) + f32(x0)
				if x_inter < x {
					winding += 1 if y0 < y1 else -1
				}
			}
		}
		if verts[i].type == VCURVE {
			x0 := i32(verts[i - 1].x)
			y0 := i32(verts[i - 1].y)
			x1 := i32(verts[i].cx)
			y1 := i32(verts[i].cy)
			x2 := i32(verts[i].x)
			y2 := i32(verts[i].y)
			ax := min(x0, min(x1, x2))
			ay := min(y0, min(y1, y2))
			by := max(y0, max(y1, y2))
			if y > f32(ay) && y < f32(by) && x > f32(ax) {
				q0: [2]f32 = {f32(x0), f32(y0)}
				q1: [2]f32 = {f32(x1), f32(y1)}
				q2: [2]f32 = {f32(x2), f32(y2)}
				hits: [2][2]f32
				if point_equal(q0, q1) || point_equal(q1, q2) {
					// Degenerate curve, treat as line
					lx0 := i32(verts[i - 1].x)
					ly0 := i32(verts[i - 1].y)
					lx1 := i32(verts[i].x)
					ly1 := i32(verts[i].y)
					if y > f32(min(ly0, ly1)) && y < f32(max(ly0, ly1)) && x > f32(min(lx0, lx1)) {
						x_inter := (y - f32(ly0)) / f32(ly1 - ly0) * f32(lx1 - lx0) + f32(lx0)
						if x_inter < x {
							winding += 1 if ly0 < ly1 else -1
						}
					}
				} else {
					num_hits := ray_intersect_bezier(orig, ray, q0, q1, q2, &hits)
					if num_hits >= 1 {
						if hits[0][0] < 0 {
							winding += -1 if hits[0][1] < 0 else 1
						}
					}
					if num_hits >= 2 {
						if hits[1][0] < 0 {
							winding += -1 if hits[1][1] < 0 else 1
						}
					}
				}
			}
		}
	}
	return winding
}

// Generate a signed distance field bitmap for a glyph
get_glyph_sdf :: proc(
	info: ^Font_Info,
	scale: f32,
	glyph: i32,
	padding: i32,
	onedge_value: u8,
	pixel_dist_scale: f32,
	width: ^i32,
	height: ^i32,
	xoff: ^i32,
	yoff: ^i32,
	allocator := context.allocator,
) -> [^]u8 {
	// Use arena allocator for intermediate allocations (optimization)
	// Stack-based backing to avoid heap allocation overhead
	arena_backing: [256 * 1024]u8
	arena: mem.Arena
	mem.arena_init(&arena, arena_backing[:])
	context.allocator = mem.arena_allocator(&arena)

	scale_x := scale
	scale_y := scale

	if scale == 0 {
		return nil
	}

	ix0, iy0, ix1, iy1: i32
	get_glyph_bitmap_box_subpixel(info, glyph, scale, scale, 0.0, 0.0, &ix0, &iy0, &ix1, &iy1)

	// If empty, return nil
	if ix0 == ix1 || iy0 == iy1 {
		return nil
	}

	ix0 -= padding
	iy0 -= padding
	ix1 += padding
	iy1 += padding

	w := ix1 - ix0
	h := iy1 - iy0

	if width != nil {
		width^ = w
	}
	if height != nil {
		height^ = h
	}
	if xoff != nil {
		xoff^ = ix0
	}
	if yoff != nil {
		yoff^ = iy0
	}

	// Invert for y-downwards bitmaps
	scale_y = -scale_y

	// Distance from singular values (in the same units as the pixel grid)
	eps: f32 : 1.0 / 1024
	eps2: f32 : eps * eps

	vertices_temp: ^Vertex = nil
	num_verts := get_glyph_shape(info, glyph, cast(^^Vertex)&vertices_temp)
	verts := ([^]Vertex)(vertices_temp)

	// Use caller's allocator for the returned bitmap
	data_ptr, _ := mem.alloc(int(w * h), allocator = allocator)
	data := ([^]u8)(data_ptr)

	// precompute uses arena (context.allocator)
	precompute_ptr, _ := mem.alloc(int(num_verts) * size_of(f32))
	precompute := ([^]f32)(precompute_ptr)

	// Precompute scaled vertex positions to avoid per-pixel multiplication
	Scaled_Vert :: struct { x, y, cx, cy: f32 }
	sv_ptr, _ := mem.alloc(int(num_verts) * size_of(Scaled_Vert))
	sv := ([^]Scaled_Vert)(sv_ptr)
	for i: i32 = 0; i < num_verts; i += 1 {
		sv[i] = {
			x  = f32(verts[i].x) * scale_x,
			y  = f32(verts[i].y) * scale_y,
			cx = f32(verts[i].cx) * scale_x,
			cy = f32(verts[i].cy) * scale_y,
		}
	}

	j := num_verts - 1
	for i: i32 = 0; i < num_verts; i += 1 {
		if verts[i].type == VLINE {
			dx := sv[j].x - sv[i].x
			dy := sv[j].y - sv[i].y
			dist := math.sqrt(dx * dx + dy * dy)
			precompute[i] = 0.0 if dist < eps else 1.0 / dist
		} else if verts[i].type == VCURVE {
			bx := sv[i].x - 2 * sv[i].cx + sv[j].x
			by := sv[i].y - 2 * sv[i].cy + sv[j].y
			len2 := bx * bx + by * by
			precompute[i] = 1.0 / len2 if len2 >= eps2 else 0.0
		} else {
			precompute[i] = 0.0
		}
		j = i
	}

	for y := iy0; y < iy1; y += 1 {
		for x := ix0; x < ix1; x += 1 {
			min_dist: f32 = 999999.0
			sx := f32(x) + 0.5
			sy := f32(y) + 0.5
			x_gspace := sx / scale_x
			y_gspace := sy / scale_y

			winding := compute_crossings_x(x_gspace, y_gspace, num_verts, verts)

			#no_bounds_check for i: i32 = 0; i < num_verts; i += 1 {
				x0 := sv[i].x
				y0 := sv[i].y

				if verts[i].type == VLINE && precompute[i] != 0.0 {
					x1 := sv[i - 1].x
					y1 := sv[i - 1].y

					dist2 := (x0 - sx) * (x0 - sx) + (y0 - sy) * (y0 - sy)
					if dist2 < min_dist * min_dist {
						min_dist = math.sqrt(dist2)
					}

					dist := abs((x1 - x0) * (y0 - sy) - (y1 - y0) * (x0 - sx)) * precompute[i]
					if dist < min_dist {
						// Check position along line
						dx := x1 - x0
						dy := y1 - y0
						px := x0 - sx
						py := y0 - sy
						t := -(px * dx + py * dy) / (dx * dx + dy * dy)
						if t >= 0.0 && t <= 1.0 {
							min_dist = dist
						}
					}
				} else if verts[i].type == VCURVE {
					x2 := sv[i - 1].x
					y2 := sv[i - 1].y
					x1 := sv[i].cx
					y1 := sv[i].cy
					box_x0 := min(min(x0, x1), x2)
					box_y0 := min(min(y0, y1), y2)
					box_x1 := max(max(x0, x1), x2)
					box_y1 := max(max(y0, y1), y2)
					// Coarse culling against bbox to avoid computing cubic unnecessarily
					if sx > box_x0 - min_dist && sx < box_x1 + min_dist && sy > box_y0 - min_dist && sy < box_y1 + min_dist {
						num: i32 = 0
						ax := x1 - x0
						ay := y1 - y0
						bx := x0 - 2 * x1 + x2
						by := y0 - 2 * y1 + y2
						mx := x0 - sx
						my := y0 - sy
						res: [3]f32 = {0, 0, 0}
						a_inv := precompute[i]
						if a_inv == 0.0 {
							// If a_inv is 0, it's 2nd degree so use quadratic formula
							qa := 3 * (ax * bx + ay * by)
							qb := 2 * (ax * ax + ay * ay) + (mx * bx + my * by)
							qc := mx * ax + my * ay
							if abs(qa) < eps2 {
								// If a is 0, it's linear
								if abs(qb) >= eps2 {
									res[num] = -qc / qb
									num += 1
								}
							} else {
								discriminant := qb * qb - 4 * qa * qc
								if discriminant >= 0 {
									root := math.sqrt(discriminant)
									res[0] = (-qb - root) / (2 * qa)
									res[1] = (-qb + root) / (2 * qa)
									num = 2
								}
							}
						} else {
							cb := 3 * (ax * bx + ay * by) * a_inv
							cc := (2 * (ax * ax + ay * ay) + (mx * bx + my * by)) * a_inv
							cd := (mx * ax + my * ay) * a_inv
							num = solve_cubic(cb, cc, cd, &res[0])
						}

						dist2 := (x0 - sx) * (x0 - sx) + (y0 - sy) * (y0 - sy)
						if dist2 < min_dist * min_dist {
							min_dist = math.sqrt(dist2)
						}

						if num >= 1 && res[0] >= 0.0 && res[0] <= 1.0 {
							t := res[0]
							it := 1.0 - t
							px := it * it * x0 + 2 * t * it * x1 + t * t * x2
							py := it * it * y0 + 2 * t * it * y1 + t * t * y2
							dist2 = (px - sx) * (px - sx) + (py - sy) * (py - sy)
							if dist2 < min_dist * min_dist {
								min_dist = math.sqrt(dist2)
							}
						}
						if num >= 2 && res[1] >= 0.0 && res[1] <= 1.0 {
							t := res[1]
							it := 1.0 - t
							px := it * it * x0 + 2 * t * it * x1 + t * t * x2
							py := it * it * y0 + 2 * t * it * y1 + t * t * y2
							dist2 = (px - sx) * (px - sx) + (py - sy) * (py - sy)
							if dist2 < min_dist * min_dist {
								min_dist = math.sqrt(dist2)
							}
						}
						if num >= 3 && res[2] >= 0.0 && res[2] <= 1.0 {
							t := res[2]
							it := 1.0 - t
							px := it * it * x0 + 2 * t * it * x1 + t * t * x2
							py := it * it * y0 + 2 * t * it * y1 + t * t * y2
							dist2 = (px - sx) * (px - sx) + (py - sy) * (py - sy)
							if dist2 < min_dist * min_dist {
								min_dist = math.sqrt(dist2)
							}
						}
					}
				}
			}

			if winding == 0 {
				min_dist = -min_dist // If outside the shape, value is negative
			}

			val := f32(onedge_value) + pixel_dist_scale * min_dist
			if val < 0 {
				val = 0
			} else if val > 255 {
				val = 255
			}
			data[(y - iy0) * w + (x - ix0)] = u8(val)
		}
	}

	// precompute and verts allocated from arena, will be freed with arena_backing
	return data
}

// Generate a signed distance field bitmap for a codepoint
get_codepoint_sdf :: proc(
	info: ^Font_Info,
	scale: f32,
	codepoint: i32,
	padding: i32,
	onedge_value: u8,
	pixel_dist_scale: f32,
	width: ^i32,
	height: ^i32,
	xoff: ^i32,
	yoff: ^i32,
	allocator := context.allocator,
) -> [^]u8 {
	return get_glyph_sdf(info, scale, find_glyph_index(info, codepoint), padding, onedge_value, pixel_dist_scale, width, height, xoff, yoff, allocator)
}

// Free a signed distance field bitmap
free_sdf :: proc(bitmap: [^]u8) {
	mem.free(bitmap)
}