package imageproc

import (
	"image"
	"math"
	"sort"

	"materialanalyzer/internal/model"
)

func ExtractCanonicalStatistics(img image.Image, visual model.VisualFeatures, texture model.TextureFeatures) model.CanonicalStatistics {
	samples := collectSamples(img, 65000)
	stats := model.CanonicalStatistics{
		MeanRGB:                 visual.AverageRGB,
		DominantOrientationDeg:  texture.PrimaryOrientationDeg,
		BrightnessHistogramBins: 16,
		GradientHistogramBins:   18,
		BrightnessDistribution:  make([]float64, 16),
		GradientHistogram:       make([]float64, 18),
	}
	if len(samples) == 0 {
		return stats
	}

	rs := make([]int, len(samples))
	gs := make([]int, len(samples))
	bs := make([]int, len(samples))
	var varianceSum float64
	for i, s := range samples {
		rs[i] = int(math.Round(s.r))
		gs[i] = int(math.Round(s.g))
		bs[i] = int(math.Round(s.b))
		dr := s.r - float64(visual.AverageRGB.R)
		dg := s.g - float64(visual.AverageRGB.G)
		db := s.b - float64(visual.AverageRGB.B)
		varianceSum += (dr*dr + dg*dg + db*db) / 3
		bin := min(len(stats.BrightnessDistribution)-1, int((s.lum/255.0)*float64(len(stats.BrightnessDistribution))))
		stats.BrightnessDistribution[bin]++
	}
	sort.Ints(rs)
	sort.Ints(gs)
	sort.Ints(bs)
	mid := len(samples) / 2
	stats.MedianRGB = model.RGB{R: rs[mid], G: gs[mid], B: bs[mid]}
	stats.ColorVariance = round4(varianceSum / float64(len(samples)) / (255.0 * 255.0))
	for i := range stats.BrightnessDistribution {
		stats.BrightnessDistribution[i] = round6(stats.BrightnessDistribution[i] / float64(len(samples)))
	}

	stats.TextureVariance, stats.GradientHistogram = gradientStatistics(img, len(stats.GradientHistogram))
	return stats
}

func gradientStatistics(img image.Image, bins int) (float64, []float64) {
	small := ResizeToMax(img, 384)
	gray, w, h := grayscale(small)
	hist := make([]float64, bins)
	if w < 3 || h < 3 || bins <= 0 {
		return 0, hist
	}

	var mags []float64
	var totalWeight float64
	for y := 1; y < h-1; y++ {
		for x := 1; x < w-1; x++ {
			gx := -int(gray[(y-1)*w+x-1]) + int(gray[(y-1)*w+x+1]) -
				2*int(gray[y*w+x-1]) + 2*int(gray[y*w+x+1]) -
				int(gray[(y+1)*w+x-1]) + int(gray[(y+1)*w+x+1])
			gy := -int(gray[(y-1)*w+x-1]) - 2*int(gray[(y-1)*w+x]) - int(gray[(y-1)*w+x+1]) +
				int(gray[(y+1)*w+x-1]) + 2*int(gray[(y+1)*w+x]) + int(gray[(y+1)*w+x+1])
			mag := math.Hypot(float64(gx), float64(gy))
			mags = append(mags, mag)
			angle := math.Atan2(float64(gy), float64(gx)) + math.Pi/2
			deg := math.Mod(angle*180/math.Pi+180, 180)
			bin := min(bins-1, int(deg/180*float64(bins)))
			hist[bin] += mag
			totalWeight += mag
		}
	}
	if len(mags) == 0 {
		return 0, hist
	}
	var mean, variance float64
	for _, mag := range mags {
		mean += mag
	}
	mean /= float64(len(mags))
	for _, mag := range mags {
		d := mag - mean
		variance += d * d
	}
	variance = variance / float64(len(mags)) / (1020.0 * 1020.0)
	if totalWeight > 0 {
		for i := range hist {
			hist[i] = round6(hist[i] / totalWeight)
		}
	}
	return round4(variance), hist
}