#include <GDIPlus.au3>
#include <Array.au3>

Opt('GUIOnEventMode', 1)

Global Const $sTitle = 'GDI+Template'
Global Const $iGuiW = 768
Global Const $iGuiH = 768
Global $bExit = False
Global $bChange = True
_GDIPlus_Startup()
Global $hGUI = GUICreate($sTitle, $iGuiW, $iGuiH + 50)
Global $hGFX = _GDIPlus_GraphicsCreateFromHWND($hGUI)
Global $hBMP = _GDIPlus_BitmapCreateFromGraphics($iGuiW, $iGuiH, $hGFX)
Global $hBUF = _GDIPlus_ImageGetGraphicsContext($hBMP)
Global $hPEN = _GDIPlus_PenCreate(0xFFE0E0E0)
_GDIPlus_GraphicsSetSmoothingMode($hBUF, 2)
Global $hSLI1 = GUICtrlCreateSlider(50, $iGuiH, $iGuiW - 50 - 100, 25), $iSLI1
Global $hLAB1 = GUICtrlCreateLabel('1000', 5, $iGuiH + 5, 50, 25, 1)
Global $hCHK1 = GUICtrlCreateCheckbox('Draw Nodes', 5, $iGuiH + 25, 100, 25)
GUICtrlCreateLabel('Knoten', $iGuiW - 100, $iGuiH + 5)
GUICtrlSetState($hCHK1, 1)
GUICtrlSetLimit($hSLI1, 1000)
GUICtrlSetData($hSLI1, 75)
GUISetOnEvent(-3, EVENT, $hGUI)
GUICtrlSetOnEvent($hCHK1, EVENT)
OnAutoItExitRegister('Xit')
GUISetState(@SW_SHOW, $hGUI)

; Return Map of Map
; [i][0] = x
; [i][1] = y
; [i][2] = t
Func _Discretize_Linear($xFunc, $iSamples)
	If $iSamples < 2 Then $iSamples = 2
	Local $Points[], $xyt[], $xy[]
	For $i = 0 To $iSamples - 1 Step 1
		$xyt[2] = $i / ($iSamples - 1)
		$xy = $xFunc($xyt[2])
		$xyt[0] = $xy[0]
		$xyt[1] = $xy[1]
		$Points[$i] = $xyt
	Next
	Return $Points
EndFunc   ;==>_Discretize_Linear

; Return Map of Map
; [i][0] = x
; [i][1] = y
; [i][2] = t
; Schlechtere Version von _Discretize_Adaptive2 was die benötigte Anzahl der Punkte angeht.
; Pro Punkt wird $xFunc aber nur genau 1x ausgewertet.
Func _Discretize_Adaptive($xFunc, $iSamples, $iIterations = 99)

	If $iSamples < 2 Then $iSamples = 2
	Local $Points = _Discretize_Linear($xFunc, Int(2 * Sqrt($iSamples))), $Map[]
	Local $iSamplesToAdd = $iSamples - UBound($Points), $OldPoints = $Map, $A, $B, $C, $p[], $j = 0, $xyt[]

	If $iSamplesToAdd <= 0 Then Return $Points

	; Da AutoIt keine Möglichkeit bietet eine performante Datenstruktur
	; wie einen RS-Tree, oder eine SortedList zu implementieren
	; wird in jedem Schritt alles kopiert...

	For $iter = 0 To $iIterations - 1 Step 1

		Local $Errors = $Map
		$Errors[0] = 0
		$Errors[UBound($Points) - 1] = 0
		For $i = 1 To UBound($Points) - 2 Step 1
			$A = $Points[$i - 1]
			$B = $Points[$i]
			$C = $Points[$i + 1]
			$p[0] = ($A[0] + $C[0]) / 2
			$p[1] = ($A[1] + $C[1]) / 2
			$Errors[$i] = Sqrt(($B[0] - $p[0]) ^ 2 + ($B[1] - $p[1]) ^ 2) ; Lokaler Error
		Next
		__ApplySmoothing($Errors)

		Local $ErrorPerEdge[UBound($Points) - 1][2] ; Arrays, weil die standard Array.au3 immernoch keine Maps supported und ich beim Standard bleiben will...
		For $i = 0 To UBound($Points) - 2 Step 1
			$ErrorPerEdge[$i][0] = ($Points[$i][2] + $Points[$i + 1][2]) / 2
			$ErrorPerEdge[$i][1] = ($Errors[$i] + $Errors[$i + 1]) / 2
		Next
		_ArraySort($ErrorPerEdge, 1, 0, 0, 1)

		; Das sind die T Werte die hinzugefügt werden. Bis zum Median der Errors.
		Local $ErrorsT[Round(UBound($ErrorPerEdge) / (2 + Sqrt($iter)))]
		For $i = 0 To UBound($ErrorsT) - 1 Step 1
			$ErrorsT[$i] = $ErrorPerEdge[$i][0]
		Next
		_ArraySort($ErrorsT)

		$j = 0
		$OldPoints = $Points
		$Points = $Map
		For $i = 0 To UBound($OldPoints) - 1 Step 1
			If $j <> -1 Then
				While $ErrorsT[$j] < $OldPoints[$i][2]
					$xyt[2] = $ErrorsT[$j]
					$xy = $xFunc($xyt[2])
					$xyt[0] = $xy[0]
					$xyt[1] = $xy[1]
					MapAppend($Points, $xyt)
					$iSamplesToAdd -= 1
					$j += 1
					If $j = UBound($ErrorsT) Or $iSamplesToAdd = 0 Then
						$j = -1
						ExitLoop
					EndIf
				WEnd
			EndIf
			MapAppend($Points, $OldPoints[$i])
		Next

		If $iSamplesToAdd = 0 Then ExitLoop

	Next

	Return $Points

EndFunc   ;==>_Discretize_Adaptive

; Return Map of Map
; [i][0] = x
; [i][1] = y
; [i][2] = t
; Bessere Version von _Discretize_Adaptive was die benötigte Anzahl der Punkte angeht.
; Dafür wird die funktion $xFunc aber deutlich häufiger ausgewertet und die Funktion ist "sehr viel langsamer"
Func _Discretize_Adaptive2($xFunc, $iSamples)

	If $iSamples < 2 Then $iSamples = 2
	Local $Points = _Discretize_Linear($xFunc, 2), $Map[], $xyt[], $xy[], $Interval[]
	Local $iSamplesToAdd = $iSamples - UBound($Points)

	Local $pOptLocal[], $pLocal[]
	Local $lX1, $lY1, $lX2, $lY2
	Local $PointsToAdd[0][0], $j = 0

	; h-Refinement
	; Starte mit einer Linie von t=0 -> t=1
	; Füge "die besten" Punkte hinzu (Punkte die den Error maximal reduzieren)
	For $iter = 0 To 199 Step 1

		$j = 0
		ReDim $PointsToAdd[UBound($Points) - 1][4] ; x, y, t, d
		For $i = 0 To UBound($Points) - 2 Step 1
			$Interval[0] = $Points[$i][2]
			$Interval[1] = $Points[$i + 1][2]
			$lX1 = $Points[$i][0]
			$lY1 = $Points[$i][1]
			$lX2 = $Points[$i + 1][0]
			$lY2 = $Points[$i + 1][1]
			$pOptLocal[2] = ($Interval[0] + $Interval[1]) / 2
			$xy = $xFunc($pOptLocal[2])
			$pOptLocal[0] = $xy[0]
			$pOptLocal[1] = $xy[1]
			$pOptLocal[3] = DistLinePoint($lX1, $lY1, $lX2, $lY2, $pOptLocal[0], $pOptLocal[1])
			For $tRange In Range('(', $Interval[0], $Interval[1], ')', 3)
				$pLocal[2] = $tRange[1]
				$xy = $xFunc($pLocal[2])
				$pLocal[0] = $xy[0]
				$pLocal[1] = $xy[1]
				$pLocal[3] = DistLinePoint($lX1, $lY1, $lX2, $lY2, $pLocal[0], $pLocal[1])
				If $pLocal[3] > $pOptLocal[3] Then
					$pOptLocal = $pLocal
				EndIf
			Next
			$PointsToAdd[$j][0] = $pOptLocal[0]
			$PointsToAdd[$j][1] = $pOptLocal[1]
			$PointsToAdd[$j][2] = $pOptLocal[2]
			$PointsToAdd[$j][3] = $pOptLocal[3]
			$j += 1
		Next

		If $iSamplesToAdd / $iSamples > 0.1 And Not $iSamplesToAdd - $iSamples = 2 Then ; Bis auf die letzten 10% werden immer mehrere Punkte hinzugefügt. Außerdem wird mit 2 Punkten gestartet.

			_ArraySort($PointsToAdd, 1, 0, 0, 3) ; Beste Punkte von oben nehmen.
			For $i = 0 To UBound($PointsToAdd) - 1 Step 1
				$PointsToAdd[$i][3] = ($i < Sqrt(UBound($PointsToAdd) - 1)) And $PointsToAdd[$i][3] > 1E-6 ; Die besten paar markieren.
			Next
			_ArraySort($PointsToAdd, 0, 0, 0, 2) ; Sortierung war ursprünglich nach t

			; Super ineffiziente Methode um einen Punkt hinzuzufügen.
			$j = 0
			$OldPoints = $Points
			$Points = $Map
			For $i = 0 To UBound($PointsToAdd) - 1 Step 1
				MapAppend($Points, $OldPoints[$i])
				If $PointsToAdd[$i][3] And $iSamplesToAdd Then
					$xyt[0] = $PointsToAdd[$i][0]
					$xyt[1] = $PointsToAdd[$i][1]
					$xyt[2] = $PointsToAdd[$i][2]
					MapAppend($Points, $xyt)
					$iSamplesToAdd -= 1
				EndIf
			Next
			MapAppend($Points, $OldPoints[UBound($OldPoints) - 1])

		Else ; Ab hier nur noch 1 Punkt hinzufügen.

			; Mark points (is it useful to add them?)
			Local $bestD = 0
			Local $bestI = 0
			For $i = 0 To UBound($PointsToAdd) - 1 Step 1
				If $PointsToAdd[$i][3] > $bestD Then
					$bestD = $PointsToAdd[$i][3]
					$bestI = $i
				EndIf
			Next

			; Super ineffiziente Methode um einen Punkt hinzuzufügen.
			$j = 0
			$OldPoints = $Points
			$Points = $Map
			For $i = 0 To UBound($PointsToAdd) - 1 Step 1
				MapAppend($Points, $OldPoints[$i])
				If $i = $bestI And $iSamplesToAdd Then
					$xyt[0] = $PointsToAdd[$i][0]
					$xyt[1] = $PointsToAdd[$i][1]
					$xyt[2] = $PointsToAdd[$i][2]
					MapAppend($Points, $xyt)
					$iSamplesToAdd -= 1
				EndIf
			Next
			MapAppend($Points, $OldPoints[UBound($OldPoints) - 1])

		EndIf

		If Not $iSamplesToAdd Then ExitLoop

	Next

	Local $A[], $B[], $C[], $bestT, $Interval2[]

	For $iter = 0 To 1 Step 1 ; 2 Iterationen reichen, wir wollen ja nicht übertreiben.

		; r-Refinement
		; Die Punkte liegen schon ganz gut (aber nicht optimal) -> Optimiere die Positionen, sodass der Error noch weiter minimiert wird
		; Im R² wird der Abstand eines Punktes zum Lotfußpunkt der Linie der Nachbarpunkte minimiert (weil das für die Anzeige am relevantesten ist, andere Fehlermaße interessieren uns nicht)
		; Um inplace zu arbeiten wird coloring verwendet. Bei einer Linie ist das supereinfach, da alle Indices Mod 2 verwendet werden können.
		; Außerdem kann die Position von t = t_min und t = t_max nicht optimiert werden weil das die Intervallgrenzen sind und die lassen wir so wie sie sind.
		For $i = 2 To UBound($Points) - 2 Step 2
			$A = $Points[$i - 1]
			$B = $Points[$i]
			$C = $Points[$i + 1]
			$Interval[0] = $A[2]
			$Interval[1] = $C[2]
			$B[3] = DistLinePoint($A[0], $A[1], $C[0], $C[1], $B[0], $B[1])
			$bestT = $B[2]
			For $tRange In Range('(', $Interval[0], $Interval[1], ')', 10)
				$pLocal[2] = $tRange[1]
				$xy = $xFunc($pLocal[2])
				$pLocal[0] = $xy[0]
				$pLocal[1] = $xy[1]
				$pLocal[3] = DistLinePoint($A[0], $A[1], $C[0], $C[1], $pLocal[0], $pLocal[1])
				If $pLocal[3] > $B[3] Then
					$B = $pLocal
					$bestT = $B[2]
				EndIf
			Next
			$Interval2[0] = $bestT - ($Interval[1] - $Interval[0]) / 12
			$Interval2[1] = $bestT + ($Interval[1] - $Interval[0]) / 12
			For $tRange In Range('(', $Interval2[0], $Interval2[1], ')', 10)
				$pLocal[2] = $tRange[1]
				$xy = $xFunc($pLocal[2])
				$pLocal[0] = $xy[0]
				$pLocal[1] = $xy[1]
				$pLocal[3] = DistLinePoint($A[0], $A[1], $C[0], $C[1], $pLocal[0], $pLocal[1])
				If $pLocal[3] > $B[3] Then
					$B = $pLocal
					$bestT = $B[2]
				EndIf
			Next
			$Points[$i] = $B
		Next

		For $i = 3 To UBound($Points) - 2 Step 2
			$A = $Points[$i - 1]
			$B = $Points[$i]
			$C = $Points[$i + 1]
			$Interval[0] = $A[2]
			$Interval[1] = $C[2]
			$B[3] = DistLinePoint($A[0], $A[1], $C[0], $C[1], $B[0], $B[1])
			$bestT = $B[2]
			For $tRange In Range('(', $Interval[0], $Interval[1], ')', 10)
				$pLocal[2] = $tRange[1]
				$xy = $xFunc($pLocal[2])
				$pLocal[0] = $xy[0]
				$pLocal[1] = $xy[1]
				$pLocal[3] = DistLinePoint($A[0], $A[1], $C[0], $C[1], $pLocal[0], $pLocal[1])
				If $pLocal[3] > $B[3] Then
					$B = $pLocal
					$bestT = $B[2]
				EndIf
			Next
			$Interval2[0] = $bestT - ($Interval[1] - $Interval[0]) / 12
			$Interval2[1] = $bestT + ($Interval[1] - $Interval[0]) / 12
			For $tRange In Range('(', $Interval2[0], $Interval2[1], ')', 10)
				$pLocal[2] = $tRange[1]
				$xy = $xFunc($pLocal[2])
				$pLocal[0] = $xy[0]
				$pLocal[1] = $xy[1]
				$pLocal[3] = DistLinePoint($A[0], $A[1], $C[0], $C[1], $pLocal[0], $pLocal[1])
				If $pLocal[3] > $B[3] Then
					$B = $pLocal
					$bestT = $B[2]
				EndIf
			Next
			$Points[$i] = $B
		Next

	Next

	Return $Points

EndFunc   ;==>_Discretize_Adaptive2

While Not $bExit And Sleep(10)

	If GUICtrlRead($hSLI1) <> $iSLI1 Then
		$iSLI1 = GUICtrlRead($hSLI1)
		GUICtrlSetData($hLAB1, $iSLI1)
		$bChange = True
	EndIf

	If $bChange Then
		_GDIPlus_GraphicsClear($hBUF)
		Local $t = TimerInit()
		Local $_ = _Discretize_Linear(F, $iSLI1)
		Draw($hBUF, Add($_, XY(0, $iGuiH / 5)), $hPEN, GUICtrlRead($hCHK1) = 1, TimerDiff($t))
		Local $t = TimerInit()
		Local $_ = _Discretize_Adaptive(F, $iSLI1)
		Draw($hBUF, Add($_, XY(0, $iGuiH / 2)), $hPEN, GUICtrlRead($hCHK1) = 1, TimerDiff($t))
		Local $t = TimerInit()
		Local $_ = _Discretize_Adaptive2(F, $iSLI1)
		Draw($hBUF, Add($_, XY(0, 4 * $iGuiH / 5)), $hPEN, GUICtrlRead($hCHK1) = 1, TimerDiff($t))
		_GDIPlus_GraphicsDrawImageRect($hGFX, $hBMP, 0, 0, $iGuiW, $iGuiH)
		$bChange = False
	EndIf

WEnd

; line (2 beliebige Punkte der Linie)
; point (Position um den Lotfußpunkt zu bestimmen)
Func DistLinePoint($lX1, $lY1, $lX2, $lY2, $pX, $pY)
	Local $vX = $lX2 - $lX1
	Local $vY = $lY2 - $lY1
	Local $pXS = $pX - $lX1
	Local $pYS = $pY - $lY1
	Local $t = ($vX * $pXS + $vY * $pYS) / ($vX ^ 2 + $vY ^ 2)
	Local $qX = $t * ($lX2 - $lX1) - $pXS
	Local $qY = $t * ($lY2 - $lY1) - $pYS
	Return ($qX ^ 2 + $qY ^ 2) ^ 0.5
EndFunc   ;==>DistLinePoint

Func __ApplySmoothing(ByRef $A)
	Local $B[]
	$B[0] = $A[0] * 0.75 + $A[1] * 0.25
	For $i = 1 To UBound($A) - 2 Step 1
		$B[$i] = $A[$i - 1] * 0.25 + $A[$i] * 0.5 + $A[$i + 1] * 0.25
	Next
	$B[UBound($A) - 1] = $A[UBound($A) - 2] * 0.75 + $A[UBound($A) - 1] * 0.25
	$A = $B
EndFunc   ;==>__ApplySmoothing

Func Draw($hGFX, $Points, $hPEN, $bDrawNodes = False, $time = 0)
	Local $xy[], $xylast[]
	For $i = 0 To UBound($Points) - 1 Step 1
		$xy = $Points[$i]
		If $bDrawNodes Then DllCall($__g_hGDIPDll, 'int', 'GdipDrawEllipse', 'handle', $hBUF, 'handle', $hPEN, 'float', $xy[0] - 2, 'float', $xy[1] - 2, 'float', 4, 'float', 4)
		If $i >= 1 Then DllCall($__g_hGDIPDll, 'int', 'GdipDrawLine', 'handle', $hBUF, 'handle', $hPEN, 'float', $xy[0], 'float', $xy[1], 'float', $xylast[0], 'float', $xylast[1])
		$xylast = $xy
	Next
	If $time And UBound($Points) > 0 Then
		_GDIPlus_GraphicsDrawString($hBUF, 'Time: ' & StringFormat('%5.2f', $time) & ' ms', 50, $Points[0][1] + 50, Default, Default, Default, 0xFFFFFFFF)

	EndIf
EndFunc   ;==>Draw

Func Range($sLeft = '[', $fMin = 0, $fMax = 1, $sRight = ']', $iSteps = 10)
	Local $A[$iSteps], $_[2], $l = $sLeft = '[' ? 0 : 1, $s = ($fMax - $fMin) / ($iSteps - 1 + ($sRight = ']' ? 0 : 1) + $l)
	For $i = 0 To $iSteps - 1 Step 1
		$_[0] = $i
		$_[1] = ($i + $l) * $s + $fMin
		$A[$i] = $_
	Next
	Return $A
EndFunc   ;==>Range

Func F($t) ; Func t€R -> [x€R, y€R]
	Local $_[]
	$_[0] = $t * $iGuiW * 0.9 + $iGuiW * 0.05
	$_[1] = 50 * ManyGauss($t * 100 - 60)
	Return $_
EndFunc   ;==>F

Func ManyGauss($x) ; Diese Funktion ist etwas stinkig. Da ist ein fieser Peak mitten in einem Tal. Linear Diskretisiert landet hier meistens maximal 'ein' Punkt.
	Return (1 * Gauss($x) + 7 * Gauss($x + 10) + 2 * Gauss($x / 5)) * Cos($x) + Gauss($x + 30) ^ 0.25 - Gauss($x + 30, 0.1)
EndFunc   ;==>ManyGauss

Func Gauss($x, $s = 1, $m = 0)
	Return Exp(-0.5 * (($x - $m) / $s) ^ 2) / ($s * 2.506628274631000502415765)
EndFunc   ;==>Gauss

Func Add($Points, $xy)
	For $i = 0 To UBound($Points) - 1 Step 1
		$Points[$i][0] += $xy[0]
		$Points[$i][1] += $xy[1]
	Next
	Return $Points
EndFunc   ;==>Add

Func XY($x, $y)
	Local $_[]
	$_[0] = $x
	$_[1] = $y
	Return $_
EndFunc   ;==>XY

Func EVENT()
	Switch @GUI_CtrlId
		Case -3
			$bExit = True
		Case $hCHK1
			$bChange = True
	EndSwitch
EndFunc   ;==>EVENT

Func Xit()
	_GDIPlus_PenDispose($hPEN)
	_GDIPlus_GraphicsDispose($hBUF)
	_GDIPlus_BitmapDispose($hBMP)
	_GDIPlus_GraphicsDispose($hGFX)
	_GDIPlus_Shutdown()
EndFunc   ;==>Xit