RUND muß es sein (Modula-2)
Polygonzüge haben manchmal einen Nachteil - es sind Polygonzüge! Das hier vorgestellte Programm nimmt dem Polygonzug seine Ecken, es interpoliert zwischen den einzelnen Punkten.
Eine Kurve ohne die von Polygonzügen typischen Ecken ist nicht ohne eine zusätzliche Berechnung möglich. Interpolation bietet sich geradezu an. um eine Kurve sichtbar zu glätten. Vielleicht werden derartige Routinen auch mal in einigen Graphikprogrammen aufgenommen!
Aber erst eine kurze Einführung in die mathematischen Notwendigkeiten. Interpolation wird benötigt, wenn man statt einer Funktionsvorschrift nur einzelne Werte kennt, beispielsweise bei einer Reihe von Meßwerten.
Mit Hilfe mathematischer Formeln werden Zwischenwerte berechnet. so. wie sie tatsächlich sein könnten. Zumindest sollte eine Interpolation so gute Werte liefern, daß sie von tatsächlichen nur minimal abweichen.
Bei unserem Programm werden ebenfalls Werte angenähert. Vorgegeben sind eine Reihe (mindestens vier) von Wertepaaren für Bildschirmkoordinaten. Durch diese Punkte wird eine Kurve gelegt. Das Beispielprogramm liefert schon eine recht anschauliche Kurve. Die Marker sollen zeigen, welche Punkte vorgegeben waren. Es sind wirklich nur zwölf!
Gerade war schon von Funktionswerten die Rede. Sie sind hier - getrennt nach x- und y-Koordinaten - die Bildschirmkoordinaten, also einerseits nacheinander die x-Koordinaten der Punkte und andererseits die y-Koordinaten. Neben diesen Werten sind noch Ableitungen nötig. Die erste Ableitung einer Funktion sagt etwas über ihr Wachstumsverhalten aus, die zweite etwas über ihr Krümmungsverhalten - und genau das interessiert uns! Wir brauchen nunmehr neben den Bildschirmkoordinaten auch das Krümmungsverhalten der Kurve in den jeweiligen Punkten.
Das ist nicht so schlimm, denn die Aufgabe der Berechnung übernimmt das Programm. Allerdings müssen zwei Werte vorgegeben werden, das Programm kann das Krümmungsverhalten am Anfangs- und Endpunkt der Kurve nicht berechnen.
Aber zwei Werte sind ja leicht anzugeben, notfalls nimmt man immer 0, was bedeutet, daß die Kurve in den beiden angesprochenen Punkten nicht gekrümmt ist.
Bevor die Kurve ausgegeben werden kann, ist die Näherung zu berechnen. Diese Aufgabe übernimmt die Routine MakeSplines. Die Ausgabe übernimmt schließlich DrawSplines. Mehr ist bei der Bedienung nicht zu beachten.
DrawSplines ist den VDI-Routinen von Megamax Modula 2 angepaßt. Daher wird die Verwendung einfacher. Wer die Graphikfunktionen des VDI nicht ausnutzen möchte, kann die Ausgabe ja auch mit LineA-Routinen abwandeln.
DrawSplines hat, wie schon angesprochen, die Aufgabe, die Kurve auszugeben. Die bereits berechneten Werten werden in die Spline-lnterpolierende (das Programm benutzt Spline-Interpolation), ein kubisches Polynom, eingesetzt und ausgewertet. Die Auswertung ist lediglich die Ausgabe der Punkte, d.h. Linien zwischen den einzelnen Punkten. Wieviele interpolierte Punkte berechnet werden sollen, ist auch in der Routine DrawSplines anzugeben. Da mit Linien gearbeitet wird, brauchen es nicht sehr viele zu sein. Je mehr Punkte berechnet werden sollen, desto länger dauert die Ausgabe auf dem Monitor.
Die beiden Routinen Square und Round werden nur intern benötigt. Da sie aber auch von allgemeinem Interesse sein könnten, sind sie auch im Definitionsmodul eingetragen.
Die Hauptaufgabe hat MakeSplines. Zuerst werden Parameter berechnet. Hat man beispielsweise eine Funktion f(x) gegeben, so existieren neben den Bildwerten f(x) auch Urbildwerte x. Die Bildschirmkoordinaten werden nur als Bildwerte angesehen, was fehlt, sind Urbildwerte.
MakeParameter liefert die Urbildwerte, einfach den Abstand der Punkte nacheinander. Danach wird in der Routine MakeSystem das Gleichungssystem aufgebaut, welches als Lösung das Krümmungsverhalten der Kurve an den einzelnen Punkten liefert, oder, um es noch etwas genauer auszudrücken, es wird die zweite Ableitung an den einzelnen Stellen berechnet. Glücklicherweise hat die Koeffizientenmatrix dieses linearen Gleichungssystems Tridiagonalgestalt, wodurch sich der Rechenaufwand verringert.
Bevor man die Koeffizienten des Interpolationspolynoms in SetSplines errechnen kann, muß noch mit SolveSystem die Lösung des Gleichungssystem ermittelt werden.
Die Koeffizienten des Interpolationspolynoms sind a,b,c und d, wobei a bei dem Term mit Grad 3 steht. In der DrawSplines-Routine ist dies deutlich zu erkennen. Wer mehr zur Theorie der Spline-Interpolation oder zu den Graphikfunktionen des VDI wissen möchte, der findet in den beiden unten angeführten Büchern sicherlich wertvolle Hinweise.
Das Listing beinhaltet keine extremen Besonderheiten von Modula-2. Eine Umsetzung in andere Programmiersprachen dürfte somit nicht schwer sein.
Literatur
[1] Numerische Mathematik, H. R. Schwarz, Teubner
[1] Atari ST Profibuch, H.-D. Jankowskn/J.F. Reschke/D. Rabich, Sybex
(***************************************************)
(* Berechnung von Splines in der Ebene zur *)
(* Glättung von Polygonzügan oder zur *)
(* Verbindung einzelner Punkte *)
(* ----------------------------------------------- *)
(* Definitionsmodul, Version 1.0, 27. April 19B8 *)
(* ----------------------------------------------- *)
(* Autor: Dietmar Rabich, Dülmen *)
(* ----------------------------------------------- *)
(* Entwickelt mit Megamax Modula 2 von Application *)
(* Systems /// Heidelberg. *)
(***************************************************)
DEFINITION NODDLE Splines;
FROM GrafBase IMPORT Point;
FRON GEMEnv IMPORT DeviceHandle;
CONST MaxPoint = 25; (* Maximalanzahl der Punkte*)
TYPE PointsArray = ARRAY [0..MaxPoint-1] OF Point; (* Punkte *)
SplinesRecord = RECORD (*Record für Splines*)
a, b, c, d,
(* Koeffizienten des Polynoms *)
XStart,XEnde : REAL (* Definitionsbereich *)
END;
SplinesArray = ARRAY [0..MaxPoint-1] OF SplinesRecord;
(* Quadrat *)
PROCEDURE Square (x : REAL) : REAL;
(* Rundung *)
PROCEDURE Round (x : REAL) : LONGINT;
(* berechnet Intarpolationapolynoa *)
PROCEDURE MakeSplines (pts : PointsArray;
NoPts : CARDINAL;
VAR splA,splB: SplinesArray;
ysa,ysb : REAL);
(* wertet Polynom aus, läuft nur unter GEM (mit VDI) *)
PROCEDURE DrawSplines (dev : DeviceHandle;
splX,splY : SplinesArray;
No : CARDINAL);
END Splines.
(* Modulende *)
Listing 1: Das Definitionsmodul
(***************************************************)
(* Berechnung von Splines in der Ebene *)
(* zur Glättung von Polygonzügan oder zur *)
(* Verbindung einzelner Punkte *)
(* ----------------------------------------------- *)
(* Implementationsmodul,Version 1.1. 24. Mai 1988 *)
(* ----------------------------------------------- *)
(* Autor: Dietmar Rabich, Dülmen *)
(* ----------------------------------------------- *)
(* Entwickelt mit Megamax Modula 2 von Application *)
(* Systems /// Heidelberg. *)
(***************************************************)
IMPLEMENTATION MODULE Splines;
(* Imports *)
FROM GrafBase IMPORT Pnt,Point;
FROM GEMEnv IMPORT DeviceHandle;
FROM MathLib0 IMPORT sqrt,power,real,entier;
FROM VDIOutputs IMPORT Line;
(* Typen *)
TYPE ParameterArray = ARRAY [0..MaxPoint-1] OF REAL; (* Parameter *)
ErgArray = ParameterArray; (* Ergebnis (y'') *)
(* Quadrat *)
PROCEDURE Square (x : REAL) : REAL;
BEGIN
RETURN power(ABS(x),2.0)
END Square;
(* Rundung *)
PROCEDURE Round (x : REAL) : LONGINT;
BEGIN
RETURN entier(x+0.5)
END Round;
(* Addition von y zu x *)
PROCEDURE IncReal (VAR x : REAL; y : REAL);
BEGIN
x:=x+y
END IncReal;
(* Routine, die Splines (Koeffizienten eines kubischen Interpolations- *)
(* polynoms) berechnet. *)
PROCEDURE MakeSplines (pts : PointsArray;
NoPts : CARDINAL;
VAR splA,splB: SplinesArray;
ysa,ysb : REAL);
TYPE GlSysArray = ErgArray;
TriGlSys = RECORD
(* u/o ; untere/obere Diagonale *)
o,u,m,l : GlSysArray
(* m : Hauptdiagonale *)
END;
(* l : Inhomogenität *)
VAR Param : ParameterArray;
SystemA,SystemB : TriGlSys;
ergA,ergB : ErgArray;
(* berechnet Parameter, damit die Koordinaten als normale Stützwerte *)
(* genutzt werden können *)
PROCEDURE MakeParameters (p : PointsArray;
VAR t : ParameterArray;
No : CARDINAL);
VAR i : CARDINAL;
BEGIN
t(0]:=0.0; (* bei 0 geht's los *)
FOR i:=1 TO No-1 DO
t[i]:=t[i-1]+ (* alter Wert + *)
sqrt(Square(real(p[i].x-p[i-1].x))+
(* Abstand der Punkte *)
Square(real(p[i].y-p[i-1].y)))
END
END MakeParameters;
(* berechnet lineares Gleichungssystem mit tridiagonaler *)
(* Koeefizientematrix *)
PROCEDURE MakeSystem (VAR system1,system2 : TriGlSys;
p : PointsArray;
t : ParameterArray;
No : CARDINAL);
VAR h : ParameterArray;
i : CARDINAL;
BEGIN
FOR i:=0 TO No-2 DO (* Schriftweiten *)
h[i]:=t[i+1]-t[i]
END;
FOR i:=0 TO No-4 DO (* Haupt- und Nebendiagonale *)
system1.a[i]:=2.0*(h[i]+h[i+1]);
system1.o[i]:=h[i+1];
system1.u[i]:=system1.o[i]
END;
system1.m[No-3]:=2.0*(h[No-3]+h[No-2]);
system2 :=system1; (* 2. System *)
FOR i:=0 TO No-3 DO (* Inhomogenität *)
system1.l[i]:=6.0*(real(p[i+1].x-p[i].x)/h[i]-real(p[i+2].x-p[i+1].x)/h[i+1]);
system2.l[i]:=6.0*(real(p[i+1].y-p[i].y)/h[i]-real(p[i+2].y-p[i+1].y)/h[i+1]);
END;
IncReal(system1.l[0], h[0] *ysa);
(* Randwarta dazu *)
IncReal(system1.l[No-3],h[No-2]*ysb);
IncReal(system2.l[0] ,h[0] *ysa);
IncReal(system2.l[No-3],h[No-2]*ysb)
END MakeSystem;
(* berechnet Lösung des Gleichungssystemes *)
(* stark vereinfacht, da Matrix Tridiagonalmatrix*)
PROCEDURE SolveSystem (system : TriGlSys;
VAR erg : ErgArray;
No : CARDINAL);
VAR i,Back : CARDINAL;
a,l,x : ErgArray;
BEGIN
m[0]:=system.m[0]; (* LR - Zerlegung *)
FOR i:=0 TO No-4 DO
l[i] :=system.u[i]/m[i]; (* L - Matrix, Nebendiagonale *)
m[i+1]:=system.m[i+1]-l[i]/system.o[i] (* R - Matrix, Hauptdiagonale *)
END;
x[0]:=system.l[0]; (* einsetzen *)
FOR i:=1 TO No-3 DO (* berechnet x aus *)
x[i]:=system.l[i]-l[i-1]*x[i-1] (* L x -system.l = 0 *)
END;
arg[No-3]:=-x[No-3]/m[No-3]; (* Lösung *)
FOR i:=0 TO No-4 DO (* berechnet erg aus *)
Back:=No-4-1; (* R erg+ x = 0 *)
erg[Back]:=-(x[Back]+system.o[Back]*erg[Back+1])/m[Back]
END
END SolveSystem;
(* berechnet die Koeffizienten des Interpolationspolynoms *)
PROCEDURE SetSplines (VAR spl1,spl2 : SplinesArray;
erg1,erg2 : ErgArray;
par : ParameterArray;
p : PointsArray;
No : CARDINAL);
VAR y : ErgArray;
i : CARDINAL;
BEGIN
(* 1. für x-Werte *)
FOR i:=0 TO No-3 DO (* 2. Ableitung *)
y[i+1]:=erg1[i]
END;
y[0] :=ysa;
y[No-1] :=ysb;
FOR i:=0 TO No-2 DO
WITH spl1[i] DO
a :=(y[i+1]-y[i])/(6.0*(par[i+1]-par[i])); (* Koeffizienten *)
b :=y[i]/2.0;
c :=real(p[i+1].x-p|i].x)/(par[i+1]-par[i])
-(par[i+1]-par[i])*(y[i+1]+2.0*y[i])/6.0;
d :=real(p[i].x);
XStart:=par[i]; (* Anfang *)
XEnde :=par[i+1] (* Ende *)
END
END;
(* 2. für y-Werte *)
FOR i:=0 TO No-3 DO (* 2. Ableitung *)
y[i+1]:=erg2[i]
END;
FOR i:=0 TO No-2 DO
WITH spl2[i] DO
a :=(y[i+1]-y[i])/(6.0*(par[i+1]-par[i])); (* Koeffizienten *)
b :=y[i]/2.0;
c :=real(p[i+1].y-p|i].y)/(par[i+1]-par[i])
-(par[i+1]-par[i])*(y[i+1]+2.0*y[i])/6.0;
d :=real(p[i].y);
XStart:=par[i]; (* Anfang *)
XEnde :=par[i+1] (* Ende *)
END
END
END SetSplines;
BEGIN
MakeParameters(pts,Param,NoPts);
MakeSystem(SystemA,SystemB,pts,Param,NoPts);
SolveSystem(SystemA,ergA,NoPts);
SolveSystem(SystemB,ergB,NoPts);
SatSplines(splA,splB,ergA,ergB,Param,pts,NoPts);
END MakeSplines;
(* wertet Interpolationspolynom aus und stellt Kurve graphisch dar *)
PROCEDURE DrawSplines (dev : DeviceHandle;
splX,splY : SplinesArray;
No : CARDINAL);
VAR i : CARDINAL;
x,Abst : REAL;
X,Y : LONGINT;
PMem,P : Point;
set : BOOLEAN;
BEGIN
set:=FALSE;
(* alle Polynome nacheinander *)
FOR i:=0 TO No-2 DO
(* alle Werte von XStart bis XEnde *)
Abst:=(splX[i].XEnde-splX[i].XStart)/20.0;
(* Feinheit der Unterteilung *)
(* ja feiner die Unterteilung, desto genauer wird die Kurve dargestellt, *)
(* aber auch desto langsamer wird die Darstellung *)
x:=splX[i].XStart;
WHILE x<=splX[i].XEnde DO
(* Interpolationapolynome *)
X:=Round(splX[i].a*power((x-splX[i].XStart),3.0) + (* für x *)
splX[i].b*Square(x-splX[i].XStart)+
splX[i].c*(x-splX[i].XStart)+
splX[i].d);
Y:=Round(splY[i].a*power((x-splY[i].XStart),3.0) + (* für y *)
splY[i].b*Square(x-splY[i].XStart)+
splY[i].c*(x-splY[i].XStart)+
splY[i].d);
P:=Pnt(SHORT(X).SHORT(Y)); (* neuer Punkt *)
IF ~set THEN
PMem:=P; (* erster Punkt als Startpunkt *)
set :=TRUE
ELSE
Line(dev,PMem,P); (* Linie zu P ziehen *)
PMem:=P
END;
IncReal(x,Abst)
END
END
END DrawSplines;
END Splines.
(* Modulende *)
Listing 2: Das Implementationsmodul
(***************************************************)
(* Berechnung von Splines in der Ebene *)
(* zur Glättung von Polygonzügen oder zur *)
(* Verbindung einzelner Punkte *)
(* ----------------------------------------------- *)
(* Demoprogramn, Version 1.0, 27. April 1988 *)
(* ----------------------------------------------- *)
(* Autor: Dietmar Rabich, Dülmen *)
(* ----------------------------------------------- *)
(* Entwickelt mit Megamax Modula 2 von Application *)
(* Systems /// Heidelberg. *)
(***************************************************)
NODULE SplinesDemo;
(* Imports *)
FRON GrafBase IMPORT Point,Pnt;
FROM AESEvents IMPORT KeyboardEvent;
FROM GEMEnv IMPORT InitGem,ExitGem,RC,
CurrGemHandle,
DeviceHandle,GemHandle;
FROM GEMGlobals IMPORT GemChar,MarkerType;
FROM VDIControls IMPORT ClearWorkstation;
FROM VDIAttributes IMPORT SetMarkerType;
FROM VDIOutputs IMPORT PolyMarker;
FROM Splines IMPORT MaxPoint,PointsArray,
SplinesRecord,
SplinesArray,
MakeSplines,DrawSplines;
VAR p : PointsArray;
(* Koordinaten der Punkte *)
n, (* Anzahl der Punkte *)
i : CARDINAL; (* Laufvariable *)
splX,splY : SplinesArray; (* Splines *)
yA,yB : REAL; (* zweite Ableitung an den Enden *)
gemHandle : GemHandle; (* Handle für GEM *)
Device : DeviceHandle; (* Workstation *)
Taste : GemChar; (* Taste *)
success : BOOLEAN; (* für Installation *)
BEGIN
InitGem(RC,Device,success);
IF success THEN (* GEM OK? *)
gemHandle:=CurrGemHandle(); (* Handle *)
yA:=0.0; (* 2. Ableitung: 0 *)
yB:=0.0; (* <-> keine Krümmung an den Enden *)
n :=12; (* 12 Punkte *)
p[0].x:=200; (* Koordinaten *)
p[0].y:=100; (* " *)
p[1].x:=100; (* " *)
p[1].y:=100; (* " *)
p[2].x:=100;
p[2].y:=200;
p[3].x:=200;
p[3).y:=200;
p[4].x:=200;
p[4].y:=300;
p[5].x:=100;
p[5].y:=300;
p[6].x:=300;
p[6].y:=100;
p[7].x:=400;
p[7].y:—100;
p[8].x:=350;
P[8].y:=100;
p[9].x:=350;
p[9].y:=200;
p[10].x:=350;
p[10].y:=250;
p[11].x:=350;
p[11].y:=300;
ClearWorkstation(Device); (* Bildschirm löschen *)
SetMarkerType(Device,crossMark);(* Kreuz als Marker *)
PolyMarker(Device,p,n);
(* Marker setzen (nicht unbedingt nötig) *)
MakeSplines(p,n,splX,splY,yA,yB);(* Splines berechnen *)
DrawSplines(Device,splX,splY,n); (* Splines ausgeben *)
KeyboardEvent(Taste); (* auf Taste warten *)
ExitGem(gemHandle) (* Auf Wiedersehen! *)
END
END SplinesDamo.
(* Demoende *)
Listing 3: Das Demoprogramm
Dietmar Rabich