ssa

1
package tg
2

3
import (
4
	"fmt"
5
	"math"
6
	"sort"
7

8
	"github.com/RB-PRO/ssa/pkg/movmean"
9
	"github.com/RB-PRO/ssa/pkg/oss"
10
)
11

12
func CalcPW(RGBs []RGB_float64, Path string) (pw []float64, Err error) {
13
	pw = make([]float64, len(RGBs))
14
	pw2 := make([]float64, len(RGBs))
15

16
	// CR
17
	for i := range RGBs {
18
		pw[i] = (RGBs[i].R*112.0 -
19
			RGBs[i].G*93.8 -
20
			RGBs[i].B*18.2) / 255.0
21

22
		// if _, err := filePW.WriteString(fmt.Sprintf("%.8f\n", pw[i])); err != nil {
23
		// 	log.Println(err)
24
		// }
25
	}
26

27
	// Вычитаем тренд
28
	pw_smooth, _ := movmean.Movmean(pw, 32)
29
	for i := range pw {
30
		pw[i] -= pw_smooth[i]
31
	}
32

33
	// Квадрат
34
	for i := range pw {
35
		pw2[i] = math.Pow(pw[i], 2)
36
	}
37

38
	// fMi := 40.0 / 60.0
39
	// cad := 30
40
	// SMO_med := cad / int(fMi)
41

42
	DEV_med := medianFilter(pw2, 30*60/40)
43
	createLineChart([]float64{}, DEV_med, Path+"DEV_med.png")
44
	// filePWmedianFilter, _ := os.Create(Path + "medianFilter.txt")
45
	// for i := range DEV_med {
46
	// 	if _, err := filePWmedianFilter.WriteString(fmt.Sprintf("%.8f\n", DEV_med[i])); err != nil {
47
	// 		log.Println(err)
48
	// 	}
49
	// }
50

51
	for i := range pw {
52
		pw[i] /= math.Sqrt(DEV_med[i])
53
	}
54
	prcMi := oss.Prctile(pw, 0.1)
55
	prcMa := oss.Prctile(pw, 99.9)
56
	prcMi = -0.0197
57
	prcMa = 0.0207
58
	fmt.Println("prcMi", prcMi, "prcMa", prcMa)
59
	createLineChart([]float64{}, pw, Path+"pw2.png")
60
	for i := range pw {
61
		if pw[i] < prcMi {
62
			pw[i] = prcMi
63
		}
64
		if pw[i] > prcMa {
65
			pw[i] = prcMa
66
		}
67
	}
68
	createLineChart([]float64{}, pw, Path+"pw.png")
69

70
	return pw, nil
71
}
72

73
func medianFilter(x []float64, n int) []float64 {
74
	// Проверка на нечетное значение n
75
	if n%2 == 0 {
76
		n++
77
	}
78

79
	// Длина входного массива
80
	length := len(x)
81

82
	// Результат фильтрации
83
	y := make([]float64, length)
84

85
	for i := 0; i < length; i++ {
86
		// Индексы для сбора значений для медианного фильтра
87
		start := i - n/2
88
		end := i + n/2
89

90
		// Гарантия, что индексы не выходят за пределы массива
91
		if start < 0 {
92
			start = 0
93
		}
94
		if end >= length {
95
			end = length - 1
96
		}
97

98
		// Извлечение значений для медианы
99
		window := x[start : end+1]
100

101
		// Сортировка окна значений и выбор медианы
102
		sortedWindow := make([]float64, len(window))
103
		copy(sortedWindow, window)
104
		sort.Float64s(sortedWindow)
105
		// Хитрый мув. При делении int(5) на int(2), получается int(2),
106
		// т.е. округление в нижнюю сторону, хотя нам нужно в старшую степень.
107
		// Поэтому из нечётного делаем чётное, а в случае получения нечётного не имеет разницы
108
		medianIndex := (len(sortedWindow) + 1) / 2
109
		// fmt.Println("medianIndex", medianIndex, "-", sortedWindow)
110
		y[i] = sortedWindow[medianIndex]
111
	}
112

113
	return y
114
}
115