sample.go

package tin

import (
	"fmt"
	"math"

	vec2d "github.com/flywave/go3d/float64/vec2"

	"github.com/flywave/go-geo"
	"github.com/flywave/go-geoid"
)

// 新增生成北京附近瓦片边界函数
func GenerateSampleTileBBox(zoom int) vec2d.Rect {
	fwLon := 118.0500 // 经度
	fwLat := 36.8167  // 纬度

	// 转换为Web墨卡托坐标
	srcProj := geo.NewProj(4326) // WGS84 经纬度
	dstProj := geo.NewProj(3857)
	pt, err := transformPoint(srcProj, dstProj, fwLon, fwLat)
	if err != nil {
		panic(err)
	}

	x, y := pt[0], pt[1]

	// 计算瓦片参数
	tileSize := (20037508.34 * 2) / math.Pow(2, float64(zoom))
	tileX := math.Floor(x / tileSize)
	tileY := math.Floor(y / tileSize)

	// 生成边界框
	return vec2d.Rect{
		Min: vec2d.T{
			tileX * tileSize,
			tileY * tileSize,
		},
		Max: vec2d.T{
			(tileX + 1) * tileSize,
			(tileY + 1) * tileSize,
		},
	}
}

// 辅助函数：坐标转换
func transformPoint(srcProj, dstProj geo.Proj, x, y float64) (vec2d.T, error) {
	points := srcProj.TransformTo(dstProj, []vec2d.T{{x, y}})
	if len(points) == 0 {
		return vec2d.T{}, fmt.Errorf("transform point failed")
	}
	return points[0], nil
}

// 创建虚拟DEM
func CreateSampleDEM(bbox vec2d.Rect, zoom int, maxElev float64) *RasterDouble {
	// 根据缩放级别动态设置分辨率
	spacing := math.Pow(2, 17-float64(zoom)) // 调整指数为17，强制16级2m、17级1m分辨率，修复层级对应问题

	cols := int(bbox.Width()/spacing) + 1
	rows := int(bbox.Height()/spacing) + 1

	data := make([]float64, cols*rows)
	for i := 0; i < rows; i++ {
		for j := 0; j < cols; j++ {
			x := bbox.Min[0] + float64(j)*spacing
			y := bbox.Min[1] + float64(i)*spacing
			// 生成波浪地形
			// 使用分形噪声生成更真实地形
			height := maxElev * (0.8*noise2D(x/1200, y/1200) +
				0.2*noise2D(x/400, y/400)) // 移除最高频噪声项，增大低频特征尺度至1200/400，进一步消除搓板效应
			data[i*cols+j] = height
		}
	}

	raster := NewRasterDoubleWithData(rows, cols, data)
	raster.SetXYPos(
		bbox.Min[0], // 左下角X坐标
		bbox.Min[1], // 左下角Y坐标（修复posy超过maxY的误差）
		spacing,
	)

	// 设置坐标转换函数（示例：Web墨卡托转WGS84）
	raster.SetTransform(func(v *Vertex) Vertex {
		// 创建原始坐标系（假设原始是Web墨卡托）
		srcProj := geo.NewProj(3857)
		// 目标坐标系（WGS84经纬度）
		dstProj := geo.NewProj(4326)

		// 执行坐标转换
		pt := srcProj.TransformTo(dstProj, []vec2d.T{{v[0], v[1]}})

		return Vertex{
			pt[0][0],
			pt[0][1],
			v[2],
		}
	})

	raster.Bounds = [4]float64{bbox.Min[0], bbox.Max[0], bbox.Min[1], bbox.Max[1]}
	return raster
}

func GenerateSampleMesh(bbox vec2d.Rect, srsCode int, zoom int) (*Mesh, error) {
	const (
		maxSurfaceElev = 50.0
	)

	// 1. 创建地理参考系
	targetProj := geo.NewProj(srsCode)
	config := &GeoConfig{
		SrcProj: targetProj,
		Datum:   geoid.HAE,
		Offset:  0,
	}

	// 2. 创建虚拟DEM并验证
	dem := CreateSampleDEM(bbox, zoom, maxSurfaceElev)
	if dem == nil || dem.Cols() == 0 || dem.Rows() == 0 {
		return nil, fmt.Errorf("DEM创建失败")
	}

	// 3. 初始化网格生成器
	zemlya := NewZemlyaMesh(config)
	if err := zemlya.LoadRaster(dem); err != nil {
		return nil, fmt.Errorf("栅格加载失败: %v", err)
	}

	// 4. 执行贪婪插入算法
	maxError := 0.5 // 最大允许误差
	zemlya.GreedyInsert(maxError)

	// 5. 转换结果网格并设置地理参考
	resultMesh := zemlya.ToMesh()

	// 6. 验证生成结果
	if resultMesh.Count() == 0 {
		return nil, fmt.Errorf("生成的网格为空")
	}
	if !resultMesh.CheckTin() {
		return nil, fmt.Errorf("生成的网格不符合TIN要求")
	}

	return resultMesh, nil
}

// 辅助哈希函数
func hash(x, y int) int {
	return (x*1619 + y*31337) & 0x7fffffff
}

// 新增二维噪声生成函数
func noise2D(x, y float64) float64 {
	// 网格单元整数部分和小数部分
	xi := int(math.Floor(x))
	yi := int(math.Floor(y))
	xf := x - float64(xi)
	yf := y - float64(yi)

	// 生成四个角点随机梯度向量
	grad := func(hash int) (float64, float64) {
		h := hash & 7
		grad := []float64{1, 1, -1, 1, 1, -1, -1, -1, 1, 0, -1, 0, 0, 1, 0, -1}[h]
		return grad, grad
	}

	// 计算四个角点贡献值
	n0, n1, n2, n3 := 0.0, 0.0, 0.0, 0.0
	t0 := 0.5 - xf*xf - yf*yf
	if t0 >= 0 {
		gx, gy := grad(hash(xi, yi))
		t0 *= t0
		n0 = t0 * t0 * (gx*xf + gy*yf)
	}

	t1 := 0.5 - (xf-1)*(xf-1) - yf*yf
	if t1 >= 0 {
		gx, gy := grad(hash(xi+1, yi))
		t1 *= t1
		n1 = t1 * t1 * (gx*(xf-1) + gy*yf)
	}

	t2 := 0.5 - xf*xf - (yf-1)*(yf-1)
	if t2 >= 0 {
		gx, gy := grad(hash(xi, yi+1))
		t2 *= t2
		n2 = t2 * t2 * (gx*xf + gy*(yf-1))
	}

	t3 := 0.5 - (xf-1)*(xf-1) - (yf-1)*(yf-1)
	if t3 >= 0 {
		gx, gy := grad(hash(xi+1, yi+1))
		t3 *= t3
		n3 = t3 * t3 * (gx*(xf-1) + gy*(yf-1))
	}

	// 归一化到[-1,1]范围
	return 40.0 * (n0 + n1 + n2 + n3)
}

type RasterAdapter struct {
	origin       *RasterDouble
	fullCoverage geo.Coverage
	tileGrid     *geo.TileGrid
	originSR     geo.Proj
}

func NewRasterAdapter(tileGrid *geo.TileGrid, coverage geo.Coverage, origin *RasterDouble) *RasterAdapter {
	return &RasterAdapter{
		origin:       origin,
		tileGrid:     tileGrid,
		fullCoverage: coverage,
		originSR:     tileGrid.Srs, // 深度克隆避免外部修改
	}
}

func intersectRect(r *vec2d.Rect, other *vec2d.Rect) *vec2d.Rect {
	minX := math.Max(r.Min[0], other.Min[0])
	minY := math.Max(r.Min[1], other.Min[1])
	maxX := math.Min(r.Max[0], other.Max[0])
	maxY := math.Min(r.Max[1], other.Max[1])

	if minX >= maxX || minY >= maxY {
		return nil
	}

	return &vec2d.Rect{
		Min: vec2d.T{minX, minY},
		Max: vec2d.T{maxX, maxY},
	}
}

func (f *RasterAdapter) generator(bbox vec2d.Rect, zoom int) (*RasterDouble, error) {
	if f.origin == nil || f.origin.Data == nil {
		return nil, fmt.Errorf("origin raster is invalid")
	}

	origin := f.origin
	// 获取分离的XY分辨率
	cellsize := origin.CellSize()
	if cellsize <= 0 {
		return nil, fmt.Errorf("invalid cellsize: %.6f", cellsize)
	}

	// 原始栅格地理范围
	originBBox := vec2d.Rect{
		Min: vec2d.T{origin.Bounds[0], origin.Bounds[2]},
		Max: vec2d.T{origin.Bounds[1], origin.Bounds[3]},
	}

	// 计算有效交集
	intersect := intersectRect(&bbox, &originBBox)
	if intersect == nil {
		return nil, fmt.Errorf("no intersection with origin raster")
	}

	// 计算列索引 (X方向)
	startCol := math.Floor((intersect.Min[0] - originBBox.Min[0]) / cellsize)
	endCol := math.Ceil((intersect.Max[0] - originBBox.Min[0]) / cellsize)

	// 计算行索引 (Y方向) - 注意Y轴从顶部开始
	startRow := math.Floor((originBBox.Max[1] - intersect.Max[1]) / math.Abs(cellsize))
	endRow := math.Ceil((originBBox.Max[1] - intersect.Min[1]) / math.Abs(cellsize))

	// 添加1个像元重叠
	overlap := 0.0
	startCol = math.Max(0, startCol)
	endCol = math.Min(float64(origin.Cols()), endCol+overlap)
	startRow = math.Max(0, startRow)
	endRow = math.Min(float64(origin.Rows()), endRow+overlap)

	// 转换为整数索引 (扩大范围确保覆盖)
	startColInt := int(math.Floor(startCol))
	endColInt := int(math.Ceil(endCol))
	startRowInt := int(math.Floor(startRow))
	endRowInt := int(math.Ceil(endRow))

	// 边界检查
	cols := origin.Cols()
	rows := origin.Rows()
	startColInt = max(0, startColInt)
	endColInt = min(cols, endColInt)
	startRowInt = max(0, startRowInt)
	endRowInt = min(rows, endRowInt)

	// 计算实际尺寸
	subCols := endColInt - startColInt
	subRows := endRowInt - startRowInt
	if subCols <= 0 || subRows <= 0 {
		return nil, fmt.Errorf("invalid subgrid size: %dx%d", subCols, subRows)
	}

	// 精确计算地理边界
	exactMinX := originBBox.Min[0] + float64(startColInt)*cellsize
	exactMaxX := originBBox.Min[0] + float64(endColInt)*cellsize
	exactMaxY := originBBox.Max[1] - float64(startRowInt)*math.Abs(cellsize)
	exactMinY := originBBox.Max[1] - float64(endRowInt)*math.Abs(cellsize)

	// === 新增：边界对齐处理 ===
	// 计算需要扩展的像元数（确保覆盖原始bbox）
	expandCols := math.Ceil((bbox.Max[0] - exactMaxX) / cellsize)
	expandRows := math.Ceil((exactMinY - bbox.Min[1]) / math.Abs(cellsize))

	if expandCols > 0 || expandRows > 0 {
		// 扩展列索引范围
		endColInt = min(origin.Cols(), endColInt+int(expandCols)+1)
		// 扩展行索引范围
		endRowInt = min(origin.Rows(), endRowInt+int(expandRows)+1)

		// 重新计算精确边界
		exactMaxX = originBBox.Min[0] + float64(endColInt)*cellsize
		exactMinY = originBBox.Max[1] - float64(endRowInt)*cellsize

		// 更新子网格尺寸
		subCols = endColInt - startColInt
		subRows = endRowInt - startRowInt
		if subCols <= 0 || subRows <= 0 {
			return nil, fmt.Errorf("expanded subgrid size invalid: %dx%d", subCols, subRows)
		}
	}

	// 高效数据复制
	totalSize := subRows * subCols
	newData := make([]float64, totalSize)
	dataSlice := origin.DataSlice()

	for r := 0; r < subRows; r++ {
		srcRow := startRowInt + r
		srcStart := srcRow*cols + startColInt
		srcEnd := srcStart + subCols
		dstStart := r * subCols

		copy(newData[dstStart:dstStart+subCols], dataSlice[srcStart:srcEnd])
	}
	const epsilon = 1e-3 // 用于处理浮点数精度的容差值

	// 创建子栅格
	subRaster := NewRasterDoubleWithData(subRows, subCols, newData)
	subRaster.SetXYPos(
		exactMinX,        // 左下角X
		exactMinY,        // 左下角Y (修复坐标计算)
		cellsize+epsilon, // 分离的分辨率
	)
	subRaster.NoData = origin.NoData
	subRaster.Bounds = [4]float64{exactMinX, exactMaxX, exactMinY, exactMaxY}

	return subRaster, nil
}

func (f *RasterAdapter) GetDEM(bbox vec2d.Rect, zoom int) (*RasterDouble, error) {
	dem, err := f.generator(bbox, zoom)
	if err != nil {
		return nil, fmt.Errorf("DEM generation failed: %w", err)
	}
	return dem, nil
}

func (f *RasterAdapter) Coverage() (geo.Coverage, error) {
	return f.fullCoverage, nil
}

// 辅助函数
func max(a, b int) int {
	if a > b {
		return a
	}
	return b
}

func min(a, b int) int {
	if a < b {
		return a
	}
	return b
}