Interpolation von werten zur kurve

[blastar]

Hallo,
ich bin ein eher seltener Schreiber da zu den meisten Fragen und Problemen schon Antworten gibt... wenn nicht hier dann in den anderen Foren.

In diesem speziellen Falle habe ich aber nichts gefunden:
Fuer ein aktuelles Projekt arbeite ich mit Messdaten welche analysiert bzw. visualisiert werden muessen. Dazu verbinde ich die vorhandenen Werte per 'Inverse Distance Weighted'... auch wenn die daraus resultierenden Daten rein spekulativ sind, die Kurve gefaellt mir gar nicht! Unten eine einfaches Beispiel, die Werte tendieren zwischen 1&2 sowie 3&4 gegen den MEAN was haessliche Beulen ergibt!

Gibt es eine Praktikable Loesung das zu harmonisieren oder einen anderen/besseren Ansatz?

Vielen Dank.

Code: Alles auswählen

Structure DATEN
  X.f
  Y.f
EndStructure
Global NewList DATEN.DATEN()

OpenWindow(0, 0, 0, 640, 490, "")
CanvasGadget(1, 0, 0, 640, 490)

AddElement(DATEN()) : DATEN()\x = 100 : DATEN()\y = 400
AddElement(DATEN()) : DATEN()\x = 300 : DATEN()\y = 400
AddElement(DATEN()) : DATEN()\x = 400 : DATEN()\y = 150
AddElement(DATEN()) : DATEN()\x = 600 : DATEN()\y = 100
  
StartDrawing(CanvasOutput(1))

  For x = 0 To 639
    a.d = 0
    b.d = 0
    ForEach DATEN()
      xd.d = (x - DATEN()\X) * (x - DATEN()\X)
      yd.d = (0 - 0) * (0 - 0)
      dist.d = xd + yd
        If dist <> 0 
          a = a + DATEN()\Y / dist
          b = b + 1 / dist
        Else
          a = a + DATEN()\Y / 1
          b = b + 1
        EndIf             
    Next
    c.d = a / b
    Circle(x, 480-c, 4, RGB(0,0,0))
  Next
  
  ForEach DATEN()
    Circle(DATEN()\X, 480-DATEN()\Y, 10, RGB(255,0,0))
  Next

StopDrawing()  

Repeat
  
  Select WaitWindowEvent()
    Case #PB_Event_CloseWindow
      End
  EndSelect    
  
ForEver

[NicTheQuick]

Das könnte dir helfen: Monotone cubic interpolation

[blastar]

Sieht ganz danach aus als ob es das ist was ich brauche ABER solche Texte bzw. Beschreibungen der Formeln habe ich noch nie verstanden und auch aus dem Schnipsel Javascript werde ich nicht schlau.

[Derren]

Ohne allzuviel Hintergrundwissen zum Thema Interpolation: Ich hätte einfach die Punkte mit Sinus Funktionen verbunden.

Sieht dann ungefähr so aus (ist gemalt, deswegen nicht ganz akkurat)



Meine Codes sind z.Z. leider ziemlich verstreut und ich finde grade nicht meinen Ansatz dazu, aber eigentlich ist es recht einfach.

Du musst nur eine Sinuskurve von #Pi/2 nach (3*#Pi)/2 zeichnen, bzw eine Cosinus-Kurve, wenn der erste Punkt unter dem zweiten liegt (oder einfach von Sinus von 3/2 Pi nach 2,5 Pi)

Wenn ich mich nicht vertan habe ist das deine Formel für abfallende Kurven:

Code: Alles auswählen

h = 1 ; "Height" - Höhendifferenz der zwei Punkte (= y1 - y2)
w = 1 ; "Width" - Horizontaler Abstand zwischen den Punkten (= x2 - x1)

h*0.5*sin(1/w* #Pi * (x+(w/2))) +(h/2) + y2

Hier Wolfram Alpha plot mit den Werten h=9, w= 10, y2=2 (d.h. y1 = 12, denn die Höhendifferenz ist 10)



Horizontale Verschiebung (das y2 am Ende ist eigtl schon die Vertikale Verschiebung):
h*0.5*sin(1/w* #Pi * (x+(w/2) + __HorizontaleVerschiebung__ )) +(h/2) + y2



Naja, ist natürlich nicht so toll, wenn du ständig fallende Werte hast, dann sieht das ja so aus....

[NicTheQuick]

Ich habe mal diesen Javascript-Code nach Purebasic übersetzt. Ist etwas hässlich geworden, aber ich habe mir auch nicht viele Gedanken über das Optimieren gemacht.
Eingefügt in deinen Ursprungscode kommt dann das hier bei raus:

Code: Alles auswählen

EnableExplicit

Structure DATEN
  X.f
  Y.f
EndStructure

Structure INTERPOLATE
	length.i
	Array xs.f(1)
	Array ys.f(1)
	Array c1s.f(1)
	Array c2s.f(1)
	Array c3s.f(1)
EndStructure
Procedure initInterpolate(List l.DATEN(), *int.INTERPOLATE)
	With *int
	\length = ListSize(l())

; 	// Rearrange xs And ys so that xs is sorted
	SortStructuredList(l(), #PB_Sort_Ascending, OffsetOf(DATEN\X), TypeOf(DATEN\X))

; 	// Get consecutive differences and slopes
	ReDim \xs.f(\length - 1)
	ReDim \ys.f(\length - 1)
	Protected Dim dys.f(\length - 1), Dim dxs.f(\length - 1), Dim ms.f(\length - 1)
; 	var dys = [], dxs = [], ms = [];
	Define i.i
	FirstElement(l())
	\xs(0) = l()\X
	\ys(0) = l()\Y
; 	For (i = 0; i < length - 1; i++) {
	For i = 0 To \length - 2
		NextElement(l())
		\xs(i + 1) = l()\X
		\ys(i + 1) = l()\Y
; 		var dx = xs[i + 1] - xs[i], dy = ys[i + 1] - ys[i];
; 		dxs.push(dx); dys.push(dy); ms.push(dy/dx);
		dxs(i) = \xs(i + 1) - \xs(i)
		dys(i) = \ys(i + 1) - \ys(i)
		ms(i) = dys(i) / dxs(i)
;	}
	Next
	
; 	// Get degree-1 coefficients
; 	var c1s = [ms[0]];
	ReDim \c1s.f(\length)
	\c1s(0) = ms(0)
; 	For (i = 0; i < dxs.length - 1; i++) {
	For i = 0 To \length - 2
; 		var m = ms[i], mNext = ms[i + 1];
		Protected m.f = ms(i)
		Protected mNext.f = ms(i + 1)
; 		If (m*mNext <= 0) {
		If (m * mNext <= 0.0)
; 			c1s.push(0);
			\c1s(i + 1) = 0
; 		} Else {
		Else
; 			var dx = dxs[i], dxNext = dxs[i + 1], common = dx + dxNext;
			Protected dx.f = dxs(i), dxNext = dxs(i + 1)
			Protected common.f = dx + dxNext
; 			c1s.push(3*common/((common + dxNext)/m + (common + dx)/mNext));
			\c1s(i + 1) = 3 * common / ((common + dxNext) / m + (common + dx) / mNext)
; 		}
		EndIf
; 	}
	Next
; 	c1s.push(ms[ms.length - 1]);
	\c1s(\length) = ms(\length - 1)
	
; 	// Get degree-2 And degree-3 coefficients
; 	var c2s = [], c3s = [];
	ReDim \c2s.f(\length - 2)
	ReDim \c3s.f(\length - 2)
; 	For (i = 0; i < c1s.length - 1; i++) {
	For i = 0 To \length - 2
; 		var c1 = c1s[i], m = ms[i], invDx = 1/dxs[i], common = c1 + c1s[i + 1] - m - m;
		common.f = \c1s(i) + \c1s(i + 1) - 2 * ms(i)
; 		c2s.push((m - c1 - common)*invDx); c3s.push(common*invDx*invDx);
		\c2s(i) = (ms(i) - \c1s(i) - common) / dxs(i)
		\c3s(i) = common / (dxs(i) * dxs(i))
; 	}
	Next
	
	EndWith
EndProcedure

Procedure.f interpolate(*int.INTERPOLATE, x.f)
	With *int
	;If (length === 0) { Return function(x) { Return 0; }; }
	If (\length = 0)
		ProcedureReturn 0.0
	EndIf
;	If (length === 1) {
	If (\length = 1)
; 		// Impl: Precomputing the result prevents problems If ys is mutated later And allows garbage collection of ys
; 		// Impl: Unary plus properly converts values To numbers
; 		var result = +ys[0];
; 		Return function(x) { Return result; };
		ProcedureReturn \ys(0)
;	}
	EndIf
; 	// Return interpolant function
; 	Return function(x) {
; 	// The rightmost point in the dataset should give an exact result
; 	var i = xs.length - 1;
	Protected i.i = \length - 1
; 	If (x == xs[i]) { Return ys[i]; }
	If (x = \xs(i))
		ProcedureReturn \ys(i)
	EndIf
;  
; 	// Search For the interval x is in, returning the corresponding y If x is one of the original xs
; 	var low = 0, mid, high = c3s.length - 1;
	Protected.i low = 0, mid, high = \length - 2
; 	While (low <= high) {
	While (low <= high)
; 		mid = Math.floor(0.5*(low + high));
		mid = (low + high) / 2
; 		var xHere = xs[mid];
		Protected xHere.f = \xs(mid)
; 		If (xHere < x) { low = mid + 1; }
		If (xHere < x)
			low = mid + 1
; 		Else If (xHere > x) { high = mid - 1; }
		ElseIf (xHere > x)
			high = mid - 1
; 		Else { Return ys[mid]; }
		Else
			ProcedureReturn \ys(mid)
		EndIf
; 	}
	Wend
; 	i = Math.max(0, high);
	i = high
	If (i < 0)
		i = 0
	EndIf
;  
;	// Interpolate
; 	var diff = x - xs[i], diffSq = diff*diff;
	Protected diff.f = x - \xs(i), diffSq.f = diff * diff
; 	Return ys[i] + c1s[i]*diff + c2s[i]*diffSq + c3s[i]*diff*diffSq;
	ProcedureReturn	\ys(i) + \c1s(i) * diff + \c2s(i) * diffSq + \c3s(i) * diff * diffSq;
; 	};
	EndWith
EndProcedure


Define NewList DATEN.DATEN()

OpenWindow(0, 0, 0, 640, 490, "")
CanvasGadget(1, 0, 0, 640, 490)

AddElement(DATEN()) : DATEN()\x = 100 : DATEN()\y = 400
AddElement(DATEN()) : DATEN()\x = 300 : DATEN()\y = 400
AddElement(DATEN()) : DATEN()\x = 400 : DATEN()\y = 150
AddElement(DATEN()) : DATEN()\x = 600 : DATEN()\y = 100

Define int.INTERPOLATE
initInterpolate(DATEN(), int)

If StartDrawing(CanvasOutput(1))
	Define x.i, y.f
	For x = 0 To 639
		y = interpolate(int, x)
		Circle(x, 480 - y, 4, RGB(0, 0, 0))
	Next
	
	ForEach DATEN()
		Circle(DATEN()\X, 480 - DATEN()\Y, 10, RGB(255, 0, 0))
	Next
	
	StopDrawing()
EndIf

Repeat
 
  Select WaitWindowEvent()
    Case #PB_Event_CloseWindow
      End
  EndSelect   
 
ForEver 

[blastar]

Derren,
danke fuer deine Muehe, hast ja einiges zusammengeschrieben (und gemalt), ist zwar nicht so ganz das was ich suche aber trotzdem interessant und irgendwann werde ich das sicher nutzen koennen!

NicTheQuick,
wow, danke fuer die schnelle Uebersetzung... auch wenn der Code schwer lesbar ist. Ich habe ihm mir mal nach und nach selber kommentiert aber so ganz steige ich da doch nicht durch... Das Grundprinzip ist mir klar (Limitierung von Anstieg und Abfall der Kurve) aber Mathe und ich stehen irgendwie auf Kriegsfuss.
Hast Du davon echt Plan davon denn leider kann ich deinen Source nciht so anpassen wie ich es mir erhofft hatte...

[NicTheQuick]

Ich habe mir die Mathematik dahinter gar nicht genau angeschaut, sondern einfach nur den Javascript-Code übersetzt und als Kommentar in meinem Code drin gelassen, sodass man sehen kann welche JS-Zeile zu welcher PB-Zeile wurde.
Eigentlich ist mein Code recht einfach einzusetzen. Ich hab ihn extra an deine vorhandene Funktion angepasst. Bloß sind aus meiner dann zwei geworden, weil ich mich so stark an den JS-Code gehalten habe. Die erste Procedure berechnet alle Hilfswerte, die man so braucht und die zweite berechnet dann die eigentliche Fitting-Kurve. Und die vorberechneten Werte werden in der Variablen int.INTERPOLATION gespeichert. Ändert man das DATEN-Array, muss man eben wieder 'initInterpolation()' aufrufen und dann kann man wieder mehrmals 'interpolate()' aufrufen.

Wie hättest du es denn gerne gehabt?

[blastar]

Urspruenglich arbeite ich mit einer XY-Map in welcher einige Daten (Messwerte) steht - das Ganze geht schon in Richtung 3D wenn man den Wert als Z betrachtet - mann sieht das es funktioniert aber die Beulen entstehen auch hier... das ganze sieht einfach nicht "rund" genug aus.
Ich hatte das ganze jetzt nur als 2D-Beispiel angegeben weil ich davon ausgegangen das es sich nur um eine einfache Erweiterung der Formel handelt welche ich dann selbst wieder umstellen koennte... scheint aber alles andere als einfach!

Code: Alles auswählen

Structure DATEN
  X.f
  Y.f
  M.f
EndStructure
Global NewList DATEN.DATEN()

OpenWindow(0, 0, 0, 600, 600, "")
CanvasGadget(1, 0, 0, 600, 600)

AddElement(DATEN()) : DATEN()\X = 49640.951 : DATEN()\Y = -70819.02 : DATEN()\M = 366.3039591969
AddElement(DATEN()) : DATEN()\X = -246.9849999999 : DATEN()\Y = -43653.42 : DATEN()\M = 368.3230506874
AddElement(DATEN()) : DATEN()\X = -25190.953 : DATEN()\Y = -125150.22 : DATEN()\M = 374.4745252922
AddElement(DATEN()) : DATEN()\X = -75078.889 : DATEN()\Y = -70819.02 : DATEN()\M = 367.472578822
AddElement(DATEN()) : DATEN()\X = -75078.889 : DATEN()\Y = 10677.78 : DATEN()\M = 330.4088411891
AddElement(DATEN()) : DATEN()\X = -246.9849999999 : DATEN()\Y = 10677.78 : DATEN()\M = 336.0524494728
AddElement(DATEN()) : DATEN()\X = 124472.855 : DATEN()\Y = 37843.38 : DATEN()\M = 323.4628594073
AddElement(DATEN()) : DATEN()\X = 49640.951 : DATEN()\Y = 119340.18 : DATEN()\M = 373.1486972655
AddElement(DATEN()) : DATEN()\X = -75078.889 : DATEN()\Y = 92174.58 : DATEN()\M = 375.1042356962
AddElement(DATEN()) : DATEN()\X = -124966.825 : DATEN()\Y = 37843.38 : DATEN()\M = 328.4530301284
AddElement(DATEN()) : DATEN()\X = 99528.887 : DATEN()\Y = -43653.42 : DATEN()\M = 364.6951314758
AddElement(DATEN()) : DATEN()\X = -25190.953 : DATEN()\Y = 92174.58 : DATEN()\M = 371.0950015584
AddElement(DATEN()) : DATEN()\X = 49640.951 : DATEN()\Y = 37843.38 : DATEN()\M = 317.899177614

Global ZMIN.f = 1000000
Global ZMAX.f = 0
ForEach DATEN()
  DATEN()\X = (DATEN()\X / 1000) + 150
  DATEN()\Y = (DATEN()\Y / -1000) + 150
  If ZMAX < DATEN()\M : ZMAX = DATEN()\M : EndIf
  If ZMIN > DATEN()\M : ZMIN = DATEN()\M : EndIf
Next

StartDrawing(CanvasOutput(1))

For y = 0 To 299
  For x = 0 To 299
    a.d = 0
    b.d = 0
    ForEach DATEN()
      xd.d = (x - DATEN()\X) * (x - DATEN()\X)
      yd.d = (y - DATEN()\Y) * (y - DATEN()\Y)
      dist.d = xd + yd
      If dist <> 0 
        a = a + DATEN()\M / dist
        b = b + 1 / dist
      Else
        a = a + DATEN()\M / 1
        b = b + 1
      EndIf             
    Next
    c.d = a / b
    d.a = (c-ZMIN) / (ZMAX-ZMIN) * 255
    d = d & %11110000
    Box(x*2,y*2,2,2,RGB(d,d,d))
  Next
Next

StopDrawing()  

Repeat
  
  Select WaitWindowEvent()
    Case #PB_Event_CloseWindow
      End
  EndSelect    
  
ForEver

[NicTheQuick]

Du könntest ja mal damit herum probieren, ob du einen schönen Skalierungswert findest, der in etwa passt.

Code: Alles auswählen

Structure DATEN
	X.f
	Y.f
	M.f
EndStructure
Global NewList DATEN.DATEN()

OpenWindow(0, 0, 0, 600, 630, "")
CanvasGadget(1, 0, 0, 600, 600)
TrackBarGadget(2, 0, 605, 550, 25, 0, 1000)
TextGadget(3, 555, 605, 40, 25, "1.0")

AddElement(DATEN()) : DATEN()\X = 49640.951 : DATEN()\Y = -70819.02 : DATEN()\M = 366.3039591969
AddElement(DATEN()) : DATEN()\X = -246.9849999999 : DATEN()\Y = -43653.42 : DATEN()\M = 368.3230506874
AddElement(DATEN()) : DATEN()\X = -25190.953 : DATEN()\Y = -125150.22 : DATEN()\M = 374.4745252922
AddElement(DATEN()) : DATEN()\X = -75078.889 : DATEN()\Y = -70819.02 : DATEN()\M = 367.472578822
AddElement(DATEN()) : DATEN()\X = -75078.889 : DATEN()\Y = 10677.78 : DATEN()\M = 330.4088411891
AddElement(DATEN()) : DATEN()\X = -246.9849999999 : DATEN()\Y = 10677.78 : DATEN()\M = 336.0524494728
AddElement(DATEN()) : DATEN()\X = 124472.855 : DATEN()\Y = 37843.38 : DATEN()\M = 323.4628594073
AddElement(DATEN()) : DATEN()\X = 49640.951 : DATEN()\Y = 119340.18 : DATEN()\M = 373.1486972655
AddElement(DATEN()) : DATEN()\X = -75078.889 : DATEN()\Y = 92174.58 : DATEN()\M = 375.1042356962
AddElement(DATEN()) : DATEN()\X = -124966.825 : DATEN()\Y = 37843.38 : DATEN()\M = 328.4530301284
AddElement(DATEN()) : DATEN()\X = 99528.887 : DATEN()\Y = -43653.42 : DATEN()\M = 364.6951314758
AddElement(DATEN()) : DATEN()\X = -25190.953 : DATEN()\Y = 92174.58 : DATEN()\M = 371.0950015584
AddElement(DATEN()) : DATEN()\X = 49640.951 : DATEN()\Y = 37843.38 : DATEN()\M = 317.899177614

Global ZMIN.f = 1000000
Global ZMAX.f = 0
ForEach DATEN()
	DATEN()\X = (DATEN()\X / 1000) + 150
	DATEN()\Y = (DATEN()\Y / -1000) + 150
	If ZMAX < DATEN()\M : ZMAX = DATEN()\M : EndIf
	If ZMIN > DATEN()\M : ZMIN = DATEN()\M : EndIf
Next

Procedure drawCanvas(scale.f = 1.0)
	Protected x.i, y.i, a.d, b.d, xd.d, yd.d, c.d, d.a
	If StartDrawing(CanvasOutput(1))
		For y = 0 To 299
			For x = 0 To 299
				a = 0
				b = 0
				ForEach DATEN()
					xd = (x - DATEN()\X) * (x - DATEN()\X)
					yd = (y - DATEN()\Y) * (y - DATEN()\Y)
					dist.d = Pow(xd + yd, scale)
					If dist <> 0
						a + DATEN()\M / dist
						b + 1 / dist
					Else
						a + DATEN()\M
						b + 1
					EndIf             
				Next
				c = a / b
				d = (c - ZMIN) / (ZMAX - ZMIN) * 255
				d = d & %11111000
				Box(x * 2, y * 2, 2, 2, RGB(d, d, d))
			Next
		Next
		StopDrawing()
	EndIf
EndProcedure

#MIN_SCALE = 0.1
#MAX_SCALE = 4.0

SetGadgetState(2, 1000.0 * (1.0 - #MIN_SCALE) / (#MAX_SCALE - #MIN_SCALE))

drawCanvas()

AddWindowTimer(0, 0, 1000)

Define scale.f = 1.0

Repeat
	
	Select WaitWindowEvent()
		Case #PB_Event_CloseWindow
			End
		Case #PB_Event_Gadget
			Select EventGadget()
				Case 2
					scale = GetGadgetState(2) * (#MAX_SCALE - #MIN_SCALE) / 1000.0 + #MIN_SCALE
					SetGadgetText(3, StrF(scale, 3))
			EndSelect
		Case #PB_Event_Timer
			If EventTimer() = 0
				drawCanvas(scale)
			EndIf
	EndSelect   
	
ForEver 

[blastar]

Leider skaliert dieser Ansatz nur die Werte, die Kurvenform selbst bleibt wie sie ist... mit den Beulen.