SB4 Interpolation

Ingenieurmathematik 1

I N G E N I E U R M A T H E M A T I K 4. Kurseinheit Autor: Prof. Dr. Edgar Jäger

Interpolation und Approximation

Die vorliegende Kurseinheit beschäftigt sich mit verschiedenen Inter-

polations- und Approximationstechniken. Zur Demonstration der

verschiedenen Interpolationstechniken werden grafische Benutzer-

oberflächen (GUI) in MATLAB erstellt.


E I NFÜHR UNG

In dieser Kurseinheit lernen Sie verschiedene Methoden zur Interpolation and Appro-

ximation kennen. Gegenstand der Interpolation sind Funktion, Kurven in Ebene und

Raum und Flächen.

Für die Bearbeitung dieses Studienbriefs sollten Sie etwa 20 Stunden einplanen.

Lernziel ��

Lernzeit ��


I NHAL T

1 Einleitung.................................................................................................................... 4

2 Lineare Interpolation................................................................................................. 6

2.1 Polynominterpolation ............................................................................................ 7

2.1.1 Berechnung des Interpolationspolynoms...................................................... 8

3 Spline-Interpolation ................................................................................................. 17

3.1 Minimaleigenschaft von Splines ......................................................................... 19

3.2 Anwendungen von Splines .................................................................................. 19

3.3 Trigonometrische Interpolation ........................................................................... 20

3.3.1 Reelle Trigonometrische Interpolation....................................................... 20

3.3.2 Komplexe Trigonometrische Interpolation ................................................ 27

3.3.3 Anwendung: Datenglättung........................................................................ 30

3.4 Bernstein-Polynome, Bézier-Darstellung von Polynomen.................................. 35

3.4.1 Bézier-Darstellung ..................................................................................... 37

3.5 Kurven................................................................................................................. 37

3.5.1 Parameterdarstellung einer Kurve .............................................................. 38

3.5.2 Parametrische kubische Splines ................................................................. 40

3.5.3 Bézier-Darstellung einer Kurve.................................................................. 42

3.6 Zweidimensionale Interpolation und Approximation.......................................... 52

3.6.1 Zweidimensionale lineare Interpolation ..................................................... 52

3.6.2 Bilineare Interpolation................................................................................ 53

3.6.3 Zweidimensionale Bikubische Splines....................................................... 54

3.6.4 Anwendung: Vergrößern digitaler Bilder................................................... 58

4 Zusammenfassung.................................................................................................... 64

5 Literatur.................................................................................................................... 65


1 E I NLE I T UNG

Eine Größe, die sich kontinuierlich ändert, ist oft nur an endlich vielen Stellen bekannt,

z. B. durch Messungen. Man möchte nun die Größe aus den zugänglichen Daten rekon-

struieren.

Interpolation kommt auch zur Anwendung, wenn ein digitales Bild vergrößert werden

soll. Ein solches Bild lässt sich durch eine Funktion B darstellen, welche einem gegebe-

nen Pixel seinen Farbwert zuordnet. Hat das Bild 100 x 100 Pixel, so ist die Funktion B

gegeben durch

B: {0,1,...,99} x {0,1,...,99} --> IR,

B(i,j) = Farbe von Pixel (i,j)

Wird nun das Bild mit Faktor 10 vergrößert, so ist diese Abbildung durch eine Abbil-

dung

IRB →× }999,...,1,0{}999,...,1,0{:~

zu ersetzen. Die naive Vorgehensweise

(1.1) ))10/(),10/((),(~

jroundiroundBjiB =

führt zu unbefriedigenden Ergebnissen:


Abbildung 1: Vergrößerung durch stückweise konstante Fortsetzung

Hier wurde eine Interpolation auf durch stückweise konstante Fortsetzung vorgenom-

men: Rundet man ab, so liefert Pixel (50,50) im Originalbild den Farbwert für die Pixel

(500+i,500+j) (i=-4,...,5, j=-4,...,5) im vergrößerten Bild.

Professionelle Foto-Bearbeitungsprogramme verwenden anstelle dieser Primitivinterpo-

lation differenzierbare Funktionen. Am häufigsten wird kubische Spline-Interpolation

eingesetzt [SCHWEIZER(2002)].


2 L I NE ARE INT E R POL AT I ON

Das einfachste Interpolationsverfahren ist die lineare Interpolation. Der Funktionswert

an einer Zwischenstelle t zwischen zwei Stützstellen 1x und 2x wird gewonnen, indem

der Wert auf der Verbindungsgeraden zwischen den gegebenen Punkten ),( 11 yx und

),( 22 yx genommen wird.

Abbildung 2: Lineare Interpolation

Hier ergibt sich also der Wert an der Zwischenstelle durch

(1.2) 12

1211 )()()(

xx

yyxtxftf

−−

−+=


2.1 POLYNOMINTERPOLATION

Ein Nachteil der linearen Interpolation sind die Knicke an den Stützstellen. Die Funkti-

on, die sich durch stückweise lineare Interpolation ergibt, ist an den Stützstellen nicht

differenzierbar:

Abbildung 3: Stückweise lineare Interpolation

Dieser Nachteil besteht bei der Polynominterpolation nicht.

SATZ (INTERPOLATIONSPLOYNOM)

Sind ),( ii yx (i=0, 1, 2,..., N) gegeben mit j)(i ≠≠ ji xx , so gibt es genau ein

Polynom p vom Grad N mit ii yxp =)( (i=0, 1, 2,..., N).


Bekanntlich ist eine Gerade (Polynom vom Grad 1) durch zwei Punkte eindeutig be-

stimmt. Obiger Satz verallgemeinert diesen Sachverhalt.

2.1.1 BERECHNUNG DES I NTERPOLATIONSPOLYNOMS

Sei NN xaxaxaaxp ++++= ...)( 2

210 .

Bei gegebenen ),( ii yx (i=0, 1, 2,..., N) ergibt sich folgendes LGS für die Koeffizien-

ten:

(2.1) iN

iNiii yxaxaxaaxp =++++= ...)( 2210 (i=0,1,...,N)

Mit der Matrix

=

−

−

−

1...

.......

.......

1...

1...

1

21

22

11

11

NNN

NN

NN

NN

xxx

xxx

xxx

V und den Vektoren

=−

0

1

.

.

a

a

a

aN

N

r,

=

Ny

y

y

y

.

.1

0

r

Mit der letzten Spalte beginnen!! Koeffizienten in umgekehrter Reihenfolge!!

ergibt sich das LGS

(2.2) yaV

rr = .

Die Matrix V heißt Vandermonde-Matrix und kann mit der MATLAB-Funktion van-

der erzeugt werden:

x-Werte der Stützstellen x=[ 1 3 4 5];

Vandermonde-Matrix

vander(x)

ans =

1 1 1 1

27 9 3 1

64 16 4 1

125 25 5 1


Die Lösung des LGS (2.2) liefert die Koeffizienten des Interpolationspolynoms.

Einfacher lässt sich das Interpolynom mit der MATLAB-Funktion polyfitpolyfitpolyfitpolyfit

bestimmen. polyfitpolyfitpolyfitpolyfit berechnet das Polynom, welches die vorgegebenen Punkte

im quadratischen Mittel am Besten annähert. Haben wir N+1 Punkte, und ist ein

Polynom vom Grad N gesucht, so erfüllt das Interpolationspolynom diese Aufgabe

exakt.

BEISPIEL 1: MATLAB-PROGRAMM „POLYNOMINTERPOLATION“ MIT

GRAFISCHER OBERFLÄCHE

Wir schreiben ein MATLAB-Programm mit grafischer Oberfläche, welches dem Benut-

zer ermöglicht, die zu interpolierenden Punkte per Mausklick einzugeben. Das Interpo-

lationspolynom, dessen Grad durch die Anzahl der vorgegebenen Punkte bestimmt ist,

wird gezeichnet.

Lösung:

Jeder Komponente eines MATLAB GUI1 (z. B. Buttons, Listboxes, Textfelder,...) ist

eine oder mehrere sog. Callback-Funktionen zugeordnet. Diese Callback-Funktionen

werden durch Aktionen des Benutzers aktiviert (z. B. Drücken eines Buttons, Maustaste

betätigen,...). Callback-Funktionen enthalten die Reaktionen auf Ereignisse [MATH-

WORKS(2006)].

1 GUI Graphical User Interface


MATLAB GUI Develop-

ment Environment starten

guide

Im Guide Quick Start

Fenster kann man zwischen

verschiedenen Vorlagen für

die die GUI-Entwicklung

wählen.

Wir wählen Blank GUI , also eine leere Vorlage


Der Bereich im Innern des

Fensters ist der sog. Layout-

Bereich. Auf diesem können

Bedienelemente angeordnet

werden. Wird die Bedien-

oberfläche gestartet, so sieht

der Benutzer diesen Bereich.

Außerhalb des Layout-

Bereichs sind Werkzeugleis-

ten angeordnet.

Wir wählen aus der Werkzeugpalette durch Anklicken das Werkzeug Achsen (Axes) und platzieren es auf dem Lay-

out-Bereich, indem wir dort den Bereich, den die Achsen einnehmen sollen, durch Ziehen festlegen.


Nun legen wir die Eigenschaften der Achsen fest. Die x-Achse soll fest auf den Bereich [0 10] eingestellt sein. Die y-

Achse soll den Bereich [-10 10] abdecken. Dazu klicken wir mit der rechten Maustaste auf den Bereich axes1 und

wählen Property Inspector. Dieser zeigt alle Eigenschaften des selektierten Bereichs an. Ferner lassen sich hier

Eigenschaften ändern.

Wir wählen die Eigenschaft XLim und stellen die Begrenzungen für x ein. Dabei wird XLimMode automatisch auf

manual umgestellt. Mit YLim verfahren wir analog.


Nun soll der Bereich, der durch die Achsen festgelegt ist, auf Mausklicks reagieren. Dazu klicken wir den Bereich

wieder mit der rechten Maustaste an und wählen View Callbacks-->ButtonDownFcn. Nun werden wir aufgefordert,

das Layout zu speichern. Gibt man einen Dateinamen ein, so wird das Layout in einer Datei mit Endung .fig gespei-

chert. Der zugehörige Programmcode befindet sich in einer Datei gleichen Namens mit Endung .m . Diese M-File

mit dem Programmcode wird nun automatisch geöffnet, und im Editor erscheint bereits die Funktion für die Behand-

lung eines Mausklicks:

% --- Executes on mouse press over axes background.

function axes1_ButtonDownFcnaxes1_ButtonDownFcnaxes1_ButtonDownFcnaxes1_ButtonDownFcn(hObject, eventdata, handles)

% hObject handle to axes1 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

Hier können wir die Anweisungen einfügen, die ausgeführt werden sollen, wenn der Benutzer auf dem entsprechen-

den Bereich eine Maustaste drückt. Zum Test geben wir hier folgende Zeile ein:

msgbox(‚hallo’);

Durch Drücken der Schaltfläche im GUI-Fenster (.fig) oder durch Drücken von im M-File Editor kann

das Programm ausgeführt werden.

Nun soll die in der Aufgabe geforderte Funktionalität implementiert werden. Dazu benötigen wir die Information, an

welcher Stelle die Maustaste gedrückt wurde. Dies ermitteln wir


mouse = get(hObject,'currentpoint');

Dabei ist hObject das Handle für das aktuelle Objekt. Dieses wird als Übergabeparameter in die Callback-

Funktion übertragen. Mit get fragen wir Eigenschaften ab. Hier handelt es sich um die Eigenschaft current-

point, welche die aktuelle Mausposition, bezogen auf die x-y-Achsen, angibt.

Mit

hold on;

plot(mouse(1,1),mouse(1,2),'ro');

wird an den angeklickten Stellen ein rotes ‚o’ gezeichnet.

Nun soll zu den angeklickten Punkten das Interpolationspolynom berechnet und gezeichnet werden. Dazu müssen

wir uns die Punkte merken. Hierfür verwenden wir globale Variable. Diese sollten beim Start des Programms initiali-

siert werden. Die Funktion xxx_OpeningFcn wird immer unmittelbar bevor die grafische Oberfläche sichtbar

wird aufgerufen. Dies ist eine geeignete Stelle für Initialisierungen:

function test1_OpeningFcOpeningFcOpeningFcOpeningFcn(hObject, eventdata, handles, varargin)

% This function has no output args, see OutputFcn.

% hObject handle to figure



% varargin command line arguments to test1 (see VARARGIN)

global points; % Vektor mit geklickten Punkten

global numberOfPoints; % Anzahl geklickter Punkte

points = zeros(2,20);

numberOfPoints = 0;

In der Mausklick-Funktion greifen wir auf diese globalen Variablen zu:

function axes1_ButtonDownFcnButtonDownFcnButtonDownFcnButtonDownFcn(hObject, eventdata, handles)




global points;

global numberOfPoints;


numberOfPoints = numberOfPoints + 1;

points(1,numberOfPoints) = mouse(1,1); % x-Koordinate

points(2,numberOfPoints) = mouse(1,2); % y-Koordinate

x = points(1,1:numberOfPoints);

y = points(2,1:numberOfPoints);


hold on;

plot(x,y,'ro');

Nun können wir das Interpolationspolynom berechnen und zeichnen:

if numberOfPoints >= 2

[p,s]=polyfit(x,y,numberOfPoints-1),

t=0:.01:10;

plot(t,polyval(p,t));

end

Unser Fenster wird nun sehr schnell mit den Interpolationspolynomen gefüllt. Daher löschen wir das vorige Interpo-

lationspolynom, ehe wir das neue zeichnen. Ferner sichern wir unser Programm noch dahingehen ab, dass bei Errei-

chen der maximalen Punktezahl von vorn begonnen wird:

function axes1_ButtonDownFcnButtonDownFcnButtonDownFcnButtonDownFcn(hObject, eventdata, handles)




global points;



numberOfPoints = numberOfPoints + 1;

points(1,numberOfPoints) = mouse(1,1); % x-Koordinate

points(2,numberOfPoints) = mouse(1,2); % y-Koordinate

x = points(1,1:numberOfPoints);

y = points(2,1:numberOfPoints);

hold on;

c = get(gca,'children'); % gca: current axes

delete( c(1:length(c)));

plot(x,y,'ro');

if numberOfPoints >= 2

[p,s]=polyfit(x,y,numberOfPoints-1),

t=0:.01:10;

plot(t,polyval(p,t));

end

if numberOfPoints == 20

msgbox('Maximalzahl Punte erreicht.');

numberOfPoints = 0;

end


Nun bauen wir noch eine Reset-Schaltfläche ein, mit der wir die Anzahl der Stützstellen auf 0 setzen können. Dazu

platzieren wir einen Pushbutton auf der Layout-Fläche, klicken auf dem Pushbutton die rechte Maustaste und wählen

Property Inspector. Eigenschaft String legt die Beschriftung der Schaltfläche fest. Hier geben wir Reset ein. Die

zugehörige Callback-Funktion ist automatisch erzeugt und kann nach Speichern editiert werden:

function pushbutton1_Callbackpushbutton1_Callbackpushbutton1_Callbackpushbutton1_Callback(hObject, eventdata, handles)

% hObject handle to pushbutton1 (see GCBO)




numberOfPoints = 0;


3 SPL I NE - I NTE R POL AT I ON

Interpolationspolynome haben die unangenehme Eigenschaft, dass sie mit zunehmen-

dem Grad und damit mit zunehmender Anzahl Stützstellen zwischen den Stützstellen i.

a. heftiger zu oszillieren neigen. Das macht Interpolationspolynome in vielen Fällen

unbrauchbar.

Man erhält i. d. R. wesentlich bessere Ergebnisse, wenn man die gegebenen Punkte

stückweise unter Verwendung von Polynomen niederen Grades interpoliert. Diese Idee

wird bei den sog. kubischen Splines verfolgt.

Spline ist die englische Bezeichnung für biegsame Holzlatten, wie sie im Schiffsbau

verwendet werden.

Bei gegebenen Punkten ),( ii yx (i=0,1,..,N) mit

Nxxx <<< ...10 setzen wir die

Spline-Funktion g auf dem Intervall [ ]1, +ii xx als Polynom dritten Grades an:

(3.1) 32)( xdxcxbaxg iiiii +++= [ ]),( 1+∈ ii xxx , i=0,1,...,N-1

Wir „kleben“ also N Polynome vom Grad drei aneinander. Dabei sollen folgende Be-

dingungen erfüllt sein:

1. iii yxg =)( (i=0,1,...,N-1)

2. 11)( ++ = iii yxg (i=0,1,...,N-1)

3. )()( 111 +++′=′

iiii xgxg (i=0,1,...,N-2)

4. )()( 111 +++″=″

iiii xgxg (i=0,1,...,N-2)

Die kubische Spline-Interpolation, die sich durch

(3.2)

[ ][ ]

[ ]

∈

∈∈

=

− ),x(x )(

...

),x(x )(

),x(x )(

)(

1-N1

211

100

NN xxg

xxg

xxg

xS


ergibt, geht also durch sämtliche vorgegebenen Punkte und ist an den Stützstellen

zweimal stetig differenzierbar.

Die Bedingungen 1. bis 4. stellen 4N-2 Bedingungen für die 4N Unbekannten Koeffi-

zienten dar. Um das System eindeutig lösbar zu machen, benötigen wir zwei weitere

Bedingungen. Man erhält somit ein LGS mit 4N Gleichungen für die 4N unbekannten

Koeffizienten.

Häufig fordert man

(3.3) 0)()( 100 =″=″− NN xgxg .

Spline-Interpolationen, die diese Zusatzbedingung erfüllen, heißen natürliche kubische

Splines.

Die MATLAB-Funktion spline verwendet anstelle von (3.3) die sog. not-a-knot-

Bedingung, bei der die äußeren drei Punkte jeweils durch ein Polynom dritten Grades

interpoliert werden:

(3.4) ])x,[x(x )()( ]),x,[x(x )()( N2-N122010 ∈=∈= −− xgxgxgxg NN

BEISPIEL: ANWENDUNG DER MATLAB-FUNKTION SPLINE

Stützstellen: Grob abgetasteter Sinus x=0:10;

y=sin(x);

Stellen, an denen Werte interpoliert werden sollen t=0:.05:10;

Splines berechnen yy=spline(x,y,t);

Ergebnis zeichnen plot(x,y,t,yy) ;


Alternativ kann man bei spline die Ableitungen an den beiden Randpunkten vorge-

ben. Ferner bietet MATLAB die Möglichkeit, mit der Funktion interp1 eine Spline-

Interpolation durchzuführen.

3.1 M INIMALEIGENSCHAFT VON SPLINES

Splines haben eine angenehme Eigenschaft, die sie für viele Anwendungen attraktiv

macht. Salopp formuliert, kommen Splines den Kurven am nächsten, die man von Hand

durch eine gegebene Punktemenge zeichnen würde. Genauer gilt für eine Spline-

Funktion folgende Eigenschaft: Ist f irgendeine hinreichend glatte Funktion, welche die

vorgegebenen Punkte interpoliert und S eine Spline-Interpolation, so gilt:

(3.5) ∫∫ ′′≤′′NN x

x

x

x

dxxfdxxS00

22)()(

Das bedeutet, dass die Spline-Funktion die kleinste mittlere Krümmung unter allen

interpolierenden Funktionen liefert. Da das obige Integral ein Maß für die Energie ist,

die erforderlich ist, um einen Stab zu biegen, kann man obige Ungleichung auch so

interpretieren, dass der Spline einen elastischen Stab mit minimaler Energie durch die

vorgegebenen Punkte bringt.

3.2 ANWENDUNGEN VON SPLINES

� Achterbahn-Konstruktion : Splines kommen bei der Achterbahn-

Konstruktion zum Einsatz. Der deutsche Achterbahn-Konstrukteur Werner

Stengel kam auf diese Idee: "Ja, in den alten Tagen wurden Achterbahnen mit

Zirkel und Lineal entworfen. Man hat einen Kreis gewählt, wo man bergab

fährt, nach dem Kreis ist man eine gerade gefahren, und nach dieser Geraden

wieder in einen Kreis. Man kann erkennen, das ist keine natürliche Bewegung,

denn wenn ich von der Geraden hier in den Kreis reinfahre, ist das eine kriti-

sche Stelle, ich werde schlagartig mit der dreifachen Belastung belastet, und

so eine Linie wie hier finden sie in der Natur gar nicht. [...]. Die Spline-

Funktionen machen so einen natürlichen Bewegungsablauf möglich, das man

z.B. beim Olympia-Looping einen Looping mit einem Glas Wasser in der Hand

fahren könnte, ohne einen Tropfen zu verschütten.“ [DEUTSCHE WEL-

LE(2005)].

� Hochgeschwindigkeitstrecken der Eisenbahn: Bei der exakten Verlegung

der Schienen kommen ebenfalls Splines zum Einsatz.


� CNC-Software enthält i. d. R. Spline-Interpolation

AUFGABE SB4-1: POLYNOM-UND SPLINEINTERPOLATION

Schreiben Sie ein MATLAB-Programm mit grafischer Oberfläche, welches dem Benutzer ermöglicht, die zu

interpolierenden Punkte per Mausklick einzugeben.

Zwei Radio Buttons sollen die Auswahl zwischen Polynominterpolation und Splineinterpolation ermögli-

chen.

3.3 TRIGONOMETRISCHE I NTERPOLATION

3.3.1 REELLE TRIGONOMETRISCHE I NTERPOLATION

Wir gehen von dem Spezialfall aus, dass 2n+1 äquidistante Stützstellen

12

2

+=

n

kxk

π (k=0,1,...,2n) mit zugehörigen Wertenky (k=0,1,...,2n)

gegeben sind. Äquidistante Stützstellen kommen in der Praxis sehr häufig vor. Wird ein

analoges Signal erfasst, so geschieht dies meist dadurch, dass das Signal mit einer be-

stimmten Rate abgetastet wird. Auf diese Weise kommen äquidistante Stützpunkte

zustande. Die Normierung auf das Intervall [0, 2π] lässt sich durch eine entsprechend

Skalierung der Zeit immer erreichen.

Gesucht ist das sog. trigonometrische Interpolationspolynom zu den vorgegebenen

Daten ),( kk yx , d. h. eine Funktion der Form

(3.6) [ ]∑=

++=n

kkk kxbkxa

axp

1

0 )sin()cos(2

)( ,

mit kk yxp =)( (k=0,1,...,2n).

��

Aufgabe

��

15 Punkte


SATZ: TRIGONOMETRISCHES INTERPOLATIONSPOLYNOM

Das trigonometrische Interpolationspolynom hat die Koeffizienten

(3.7a) ∑=+

=n

jjjk kxy

na

2

0

)cos(12

2 (k=0,1,...,n)

(3.7b) ∑=+

=n

jjjk kxy

nb

2

0

)sin(12

2 (k=1,2,...,n)

Der Algorithmus von Goertzel [VERFÜRTH(2005)] gestattet eine effiziente Berech-

nung der Koeffizienten und eine effiziente Auswertung des trigonometrischen Interpola-

tionspolynoms.

ALGORITHMUS VON GOERTZEL

Berechne

∑=+

=n

jjy

na

2

00 12

2

Für k=1,...,n:

)sin(,2),cos( kk xscccxc ===

Mit 02212 == ++ nn uu berechne für j=2n, 2n-1,...,1:

21 ++ −⋅+= jjjj uuccyu

Damit:

( )

1

210

12

212

2

usn

b

uucyn

a

k

k

⋅⋅+

=

−⋅++

= (k=1,...,n)


Auswertung des trigonometrischen Polynoms an einer Stelle *x :

Setze

)sin(,2),cos(,0 **21 xscccxcuu nn ===== ++

Berechne für j=n,n-1,...,1:

21 ++ −⋅+= jjjj uuccau

Dann gilt: 1210*

2

1)( usuucaxp ⋅+−⋅+=

MATLAB-Code für trigonometrische Interpolation

Berechnung der Koeffizienten: Datei trigoInt.m

function koeff = trigoInttrigoInttrigoInttrigoInt(y)

%

% Trigonometrische Interpolation an den Stützstellen

% xk=2*pi*k/(2n+1), k=0,...,2n

%

% Eingabeparameter: y : Zeilenvektor der Werte an den Stellen xk

% Länge von y: 2n+1

%

% Rückgabewerte: koeff: Koeffizienten des trigonometrischen

% Interpolationspolynoms

% 1. Zeile: a0,a1,...,an

% 2. Zeile: 0, b1,...,bn

%

nn=size(y);

m=nn(2); % m = 2n+1; m-1=2n

n = (m-1)/2;

h=2*pi/m; % Abstand der Stützstellen

x=zeros(1,m);


for i=2:m

x(i)=h*(i-1); % Stützstellen

end

hold off;

plot(x,y,'ro');

hold on;

z=x(1,2:n+1);

c=cos(z);

cc=2*c;

s=sin(z);

u = zeros(n,m+1);

% a0 = koeff(1,1)

% a1 = koeff(1,2),...

% b1 = koeff(2,2), b2 =koeff(2,3),...

koeff=zeros(2,n+1);

koeff(1,1)=0; % a0 = koeff(1,1)

for j=1:m

koeff(1,1)=koeff(1,1)+y(j);

end

koeff(1,1)=koeff(1,1)*2/m; % a0

for k=1:n

for j=m-1:-1:1

u(k,j)=y(j+1)+cc(k)*u(k,j+1)-u(k,j+2);

end

koeff(1,k+1)=(y(1)+c(k)*u(k,1)-u(k,2))*2/m; %ak

koeff(2,k+1)=s(k)*u(k,1)*2/m; % bk

end

Auswertung der trigonometrischen Interpolation: Datei trigoIntValue.m

function values = trigoIntValuetrigoIntValuetrigoIntValuetrigoIntValue(x,y) % % Berechnet den Wert der trigonometrischen Interpolation % an der Stelle x % % Eingabeparameter: x: Stellen, an denen interpoliert werden soll % (Zeilenvektor) % y: Zeilenvektor der Werte an den Stellen % xk=2*pi*k/(2n+1), k=0,...,2n % Länge von y: 2n+1


% % Rückgabewerte: Werte an den Stellen x % koeff = trigoInt(y); nn=size(y); m=nn(2); % m = 2n+1 n = (m-1)/2; ww=size(x); w=ww(2); u=zeros(w,n+2); c=cos(x); cc=2*c; s=sin(x); values=zeros(ww); r=zeros(ww); for i=1:w for j=n:-1:1 u(i,j)=koeff(2,j+1)+cc(i)*u(i,j+1)-u(i,j+2); end r(i)=s(i)*u(i,1); for j=n:-1:1 u(i,j)=koeff(1,j+1)+cc(i)*u(i,j+1)-u(i,j+2); end values(i) = 0.5*koeff(1,1)+c(i)*u(i,1)-u(i,2)+r(i); end

ANWENDUNGSBEISPIELE

>> h=2*pi/5;

>> xx=0:h:2*pi-h;

>> y=sin(xx)+cos(xx);

>> x=0:.01:2*pi;

>> plot(x,trigoIntValue(x,y))


>> y=rand(1,5);

>> x=0:.01:2*pi;

>>plot(x,trigoIntValue(x,y));

>> y=rand(1,21);

>> x=0:.01:2*pi;

>>plot(x,trigoIntValue(x,y));


Ist die Anzahl der Stützstellen eine Zweierpotenz, so lässt sich die trigonometrische

Interpolation besonders effizient mit der sog. schnellen Fourier-Transformation (FFT,

Fast Fourier Transform) berechnen. Diese wird mit der MATLAB-Funktion fft

durchgeführt.

AUFGABE SB4-2: HOCH-UND TIEFPASSFILTER

Schreiben Sie ein MATLAB-Programm mit grafischer Oberfläche, welches zunächst die

Funktion

f(x) = x (0 ≤ x ≤ 0.5), f(x) = x - 0.5 (0.5< x ≤ 1)

an 101 Stützstellen trigonometrisch interpoliert.

Der Benutzer soll nun die minimale Frequenz und die maximale Frequenz eingeben

können. Entsprechend der Eingabe werden die Koeffizienten (3.7) des trigonometrischen

Polynoms 0 gesetzt und das trigonometrische Interpolationspolynom (3.6) geplottet.

Dazu ist insbesondere eine der Funtionen trigoInt() oder trigoIntValue() zu

modifizieren.

��

Aufgabe

��

15 Punkte


3.3.2 K OMPLEXE TRIGONOMETRISCHE I NTERPOLATION

Die MATLAB-Funktion fft geht von der komplexen Formulierung des

trigonometrischen Interpolationsproblems aus:

KOMPLEXE FORMULIERUNG DER TRIGONOMETRISCHEN INTERPO -

LATION

Gegeben sind n äquidistante Stützpunkte ( )jj yx , (j=0,1,...,n-1) mitn

Tjx j = . Die

jy können komplexe Zahlen sein.

Gesucht ist eine Funktion der Form

kxin

kkeY

nxy

α−−

=∑=

1

0

1)( ,

welche an den Stützstellen mit den jy übereinstimmt.

Nimmt man an, dass die n gegebenen Stützpunkte Teil einer periodischen Funktion mit

der Periode T sind, so ergibt sich aus der Forderung, dass auch y(x) T-periodisch sein

soll,T

πα 2= . Die Fourier-Koeffizienten sind gegeben durch

∑−

==

1

0

2n

j

n

kji

jk eyYπ

(k=0,1,...,n-1)

Die Berechnung der kY aus den jy nennt man Fourier-Transformation (FT). Die

Auswertung von

n

kjin

kkj eY

nxy

π21

0

1)(

−−

=∑=

bei gegebenen kY heißt inverse Fourier-Transformation (IFT ).


3.3.2.1 Spezialfall: Reelle Stützwerte

Hat man reelle Stützwerte jy gegeben, so ergibt sich

∑∑−

=

−

=

+==1

0

1

0

2 2sin

2cos

n

jj

n

j

n

kji

jk n

kji

n

kjyeyY

πππ

,

und für die konjugiert komplexen Fourier-Koeffizienten

∑−

=

−=1

0

2sin

2cos

n

jjk n

kji

n

kjyY

ππ

Damit folgt:

k

n

jj

n

jjkn

Yn

kji

n

kjy

n

jkni

n

jknyY

=

−−

−=

=

−−−=

∑

∑−

=

−

=−

1

0

1

0

22sin

22cos

)(2sin

)(2cos

ππππ

ππ

Es gilt also die Symmetriebedingung

(S) kkn YY =− (k=1,...,n-1)

Für n=5 gilt bedeutet dies z. B.

413241 ,, YYYYYY === ,

und für n=6 ergibt sich

334251 ,, YYYYYY === ,

also ist in diesem Fall 3Y reell.

Sind umgekehrt die Stützwerte jy reell und gilt (S), so ist die IFT reell:


Mit ni

eπ

ω2

= gilt:

( )z

n

ki

n

keeez kn

ki

ikn

ni

kn =−+−===== −−−− )2

sin()2

cos(2

2)

2 ππωωπ

ππ

.

Fall1: n gerade, also n = 2m. Dann ist mY reell, und es folgt

[ ]

( ) ( ) IRYn

jkYYi

n

jkYYY

n

Yn

jki

n

jkY

n

jki

n

jkYY

n

YYYYn

YYYYYYYYn

YYYYYYYYn

YYYYY

YYYn

Yn

eYn

xy

jm

m

kkkkk

mjm

m

kk

jk

k

mjm

m

k

jkk

jkk

mjm

jmm

jmm

jjjj

mjm

jmm

jmm

jn

jjn

j

mjm

jmm

jmm

njn

j

njn

jn

k

kjk

n

kjin

kkj

∈−+

+−+++=

+

++

−+=

+++=

++++++++=

++++++++=

+++++

+++===

−−

=

−−

=

−

−−

=

−

−−−

−−−

−−

−−+

−−−−

−−

−

−+−+

−−−

−−−

−

−−−

−−

=

−−−

=

∑

∑

∑

∑∑

])1()2

sin()2

cos([1

])2

sin()2

cos()2

sin()2

cos([1

][1

]...[1

]...[1

]...

[111

)(

1

10

1

10

1

10

)1(1

)1(1

22

22110

)1(1

)1(1

22

22110

)1(1

)1(1

)2(2

22

)1(110

1

0

21

0

ππ

ωππππ

ωωω

ωωωωωωω

ωωωωωωω

ωωωωω

ωωωπ

Fall 2: n ungerade, also n = 2m+1. Dann folgt


[ ]

( ) ( ) IRn

jkYYi

n

jkYYY

n

n

jki

n

jkY

n

jki

n

jkYY

n

YYYn

YYYYYY

YYYn

Yn

xy

m

kkkkk

m

kk

jk

k

m

k

jkk

jkk

jmm

mjm

jmm

jmm

njn

j

njn

jm

k

kjkj

∈

+−+++=

++

−+=

++=

++++++

+++==

∑

∑

∑

∑

=

=

−

=

−

+−+

−+−+

−−−

−−−

−

−−−

−

=

−

])2

sin()2

cos([1

])2

sin()2

cos()2

sin()2

cos([1

][1

]...

[11

)(

10

10

10

)1(1

)2(2

)1(1

)2(2

22

)1(110

2

0

ππ

ππππ

ωω

ωωωωωω

ωωω

3.3.3 ANWENDUNG : DATENGLÄTTUNG

Als Beispiel betrachten wir folgende Funktion im Intervall [0, 4]:

(#)

++= −−−−−)10sin(

3

132)(

22 )3()1(22 xeeexf xxx

[SORMANN(2006)].

function f=testfn(x) f=exp(-x/2).*(2*exp(-2*(x-1).^2)+3*exp(-2*(x-3).^2)+sin(10*x)/3);


Abbildung 4: Testfunktion f(x) (#)

Wir überlagern nun obige Funktion mit zufälligen Störungen:

>> h=4/64;

>> xx=0:h:4-h;

>> y=testfn(xx);

>> z=y+.1-0.2*rand(1,64)


Abbildung 5: Funktion (#) mit Störungen

FFT für ungestörtes Signal und für gestörtes Signal, jeweils an 64 äquidistanten Stützstellen abgetastet:

>> yClean = testfn(xx); % Werte ungestört

>> uClean = fft(yClean); % FFT für ungestörte Werte

>> plot(abs(uClean));

>> uDirty = fft(z); % FFT für gestörte Werte

>> plot(abs(uDirty),'r');


Abbildung 6: Beträge der komplexen Fourier-Koeffizienten

Wie man sieht, wirken sich die Zufallsstörungen so aus, dass im mittleren Frequenzbereich von Null verschiedene

Koeffizienten auftreten. Es liegt nahe, das Signal durch Unterdrücken der mittleren Frequenzen zu glätten:

>>k0=12;

>>for i=k0+1:64-k0+1

>> uDirty(i)=0

>>end

>> yCleaned = ifft(uDirty); % inverse FFT

>> plot(real(yCleaned),’r’); % Realteil der inversen FT

>> hold on;

>> plot(y,’r’); % zum Vergleich: Original-Funktion


Abbildung 7: Original-Funktion, gestörte Werte und Glättung


3.4 BERNSTEIN-POLYNOME , BÉZIER -DARSTELLUNG VON POLYNOMEN

Nach dem binomischen Lehrsatz gilt:

( )( ) ∑∑==

− =−

=+−=

n

kkn

n

k

kknn tBttk

ntt

00

)()1(11

mit den Bernstein-Polynomen vom Grad n bezüglich des Intervalls [0, 1]

kknkn tt

k

ntB −−

= )1()(

Bernstein-Polynome haben folgende Eigenschaften:

1. )(tBkn hat in [0, 1] genau ein Maximum, und zwar bei n

kt =max .

2. t=0 ist k-fache Nullstelle von )(tBkn

3. t=1 ist (n-k)-fache Nullstelle von )(tBkn

4. Die Bernstein-Polynome genügen der Rekursionsformel

)()(

)()1()(

)()1()()(

1,1

1,00

1,1,1

ttBtB

tBttB

tBtttBtB

nnnn

nn

nknkkn

−−

−

−−−

=−=

−+=


BEISPIEL: BERNSTEIN-POLYNOME FÜR N=4

( ) ( ) ( ) ( )

)()1()()(

)()1()()(

)()1()()(

),()1()()(

)()1()()(

,)1()()1()()(

)(,)(,)(,)(

1)(,1)(,1)(,1)(,1)(

332334

231324

120223

130314

120213

2110112

444

333

22211

001

012

023

034

04

tBtttBtB

tBtttBtB

tBtttBtB

tBtttBtB

tBtttBtB

ttttBtttBtB

ttBttBttBttB

tBttBttBttBttB

−+=−−=−+=−+=−+=

+−=−+=

====

=−=−=−=−=

Abbildung 8: Berstein-Polynome 4. Grades


3.4.1 BÉZIER -DARSTELLUNG

Jedes Polynom p vom Grad n lässt sich als Linearkombination der Bernstein-Polynome

darstellen:

(B) )()(0

tBtp kn

n

kk∑

=

= β

(B) heißt Bézier2-Darstellung von p. Die kβ heißen Bézier-Koeffizienten und die

Punkte (i/k, kβ ) Bézier-Punkte.

3.5 K URVEN

In CAD-Anwendungen hat man häufig das Problem, Kurven darzustellen, die durch

vorgegebene Punkte gehen müssen.

Industrieroboter müssen häufig möglichst genau Bahnen (Trajektorien) verfolgen,

wobei zusätzlich eine konstante Bahngeschwindigkeit gefordert sein kann. Dies ist z.

B. beim Bahnschweißen, Entgraten und beim Kleberauftrag der Fall.

Die Bewegung kann von Punkt zu Punkt erfolgen, d. h. es ist eine Reihe von Raum-

punkten

( )iii zyx ,, (i=1,2,...,k)

zu vorgegebenen Zeitpunkten kttt ,...,, 21 mit vorgegebenen Geschwindigkeiten

kvvv ,...,, 21 zu durchlaufen. Die Bewegung soll dabei ruckfrei sein, d. h. die zweite

Ableitung der Interpolationsfunktion (Beschleunigung) soll stetig sein.

2 Bézier, Pierre: Französischer Mathematiker (1910-1999). Beschrieb 1960 die nach

ihm benannte Bézier-Kurve, die von seinem damaligen Arbeitgeber Renault zur Gestal-

tung von Karosserieformen eingesetzt wurde.


3.5.1 PARAMETERDARSTELLUNG EINER K URVE

Eine Parameterdarstellung einer Kurve in der Ebene bzw. im Raum hat die Form

=

)(

)()(

ty

txtr

r bzw.

=)(

)(

)(

)(

tz

ty

tx

trr

BEISPIEL: PARAMETERDARSTELLUNGEN VON KURVEN

=

)sin(5

)cos(5

)(

)(

t

t

ty

tx

Kreis um (0,0) mit Radi-

us 5

>>t=0:.01:10;

>> plot(5*cos(t),5*sin(t))

=

−

−

)sin(5

)cos(5

)(

)(1.0

1.0

te

te

ty

txt

t

Logarithmische Spirale

>> t=0:.01:40;

>> plot(5*exp(-0.1*t).*cos(t),5*exp(-0.1*t).*sin(t));


−+

−=

3

2

2

1

0

1

)(

)(

)(

t

tz

ty

tx

Gerade

>> t=0:.01:20;

>> plot3(1+2*t,-2*t,-1+3*t);


=

t

t

t

tz

ty

tx

)sin(

)cos(

)(

)(

)(

zylindrische Spirale (He-

lix, Schraubenlinie)

>> t=0:.01:40;

>> plot3(cos(t),sin(t),t);

3.5.2 PARAMETRISCHE KUBISCHE SPLINES

Die oben diskutierten Spline-Funktionen können eingesetzt werden, um eine glatte

Kurve durch eine Reihe von gegebenen Punkten zu legen, wobei auch die Zeitpunkte,

zu denen die Punkte passiert werden sollen, gegeben sein können. Dabei werden einfach

Splines für jede Komponente berechnet.

BEISPIEL: PUNKTE AUF DER ZYLINDRISCHEN SPIRALE VORG EGEBEN

Gesucht ist eine parametrische Spline-Kurve

)(

)(

)(

tz

ty

tx

mit

=

i

i

i

i

i

i

t

t

t

tz

ty

tx

)sin(

)cos(

)(

)(

)(

>> % vorgegebene Punkte

>>t=0:5:40 ;

>> x=cos(t);

>> y=sin(t);


>> z=t;

>>tt=0:.01:40; % an diesen Stellen sollen die Kurvenpunkte berechnet werden

>> xx = spline(t,x,tt);

>> yy=spline(t,y,tt);

>> zz=spline(t,z,tt);

>> plot3(xx,yy,zz,'r');

>> hold on;

>> plot3(x,y,z,'*');

>>


AUFGABE SB4-3: SPLINE-INTERPOLATIONSKURVE

In der Datei points.matpoints.matpoints.matpoints.mat sind 26 Punkte ( )ii yx , (i=1,2,...,26) binär gespeichert. Mit

load points.matload points.matload points.matload points.mat werden die Variablen x und y mit den entsprechenden Werten

geladen.

In der Ebene soll eine Bahn abgefahren werden, bei der die Punkte in der gegebenen

Reihenfolge durchlaufen werden. Dabei soll der i-te Punkt nach

2

25

110

−= it i Se-

kunden durchlaufen werden (i=1,…,26), so dass der letzte Punkt nach 10 Sekunden er-

reicht ist. Die Bahn zwischen den Punkten soll durch eine Spline-Kurve gegeben sein.

An welcher Stelle befindet man sich nach 5 Sekunden?

3.5.3 BÉZIER -DARSTELLUNG EINER K URVE

Sind n+1 Kontrollpunkte3 2IRbi ∈r

bzw. 3IRbi ∈r

gegeben, so heißt

∑=

=n

kknk tBbtP

0

)()(rr

die durch die Punkte ibr

erzeugte Bézier-Kurve. Die )(tBkn sind dabei die Bernstein-

Polynome n-ten Grades.

Es gilt:

n

n

kknk

n

kknk bBbPbBbP

rrrrrr==== ∑∑

== 00

0

)1()1( ,)0()0(

Die Bézier-Kurve interpoliert also die Endpunkte. Ferner zeigt der Tangentialvektor in

den Endpunkten in Richtung des Verbindungsvektors der ersten bzw. letzten beiden

Kontrollpunkte.

3 Eigentlich: Ortsvektoren von Kontrollpunkten

��

Aufgabe

��

10 Punkte


Im Fall von 2 Kontrollpunkten ergibt die Bézier-Kurve die Verbindungsstrecke der

beiden Punkte. Bei drei Kontrollpunkten liegt die Bézier-Kurve im Innern des von den

drei Punkten ausgespannten Dreiecks. Diese Eigenschaft lässt sich auf beliebig viele

Kontrollpunkte verallgemeinern: Die Bézier-Kurve bewegt sich in der sog. konvexen

Hülle der Kontrollpunktmenge.

Wird ein Kontrollpunkt geändert, so ändert sich in der Bézier-Kurve nur der zugehörige

Summand.

3.5.3.1 Der Algorithmus von de Casteljau

Aus der Rekursion für die Bernstein-Polynome lässt sich ein Algorithmus zu r Berech-

nung von Bézier-Kurven ableiten (z. B. [BABOVSKY, NEUNDORF(2006)]:


ALGORITHMUS VON DE CASTELJAU

Gegeben seien die n+1 Kontrollpunkte 2IRbi ∈r

bzw. 3IRbi ∈r

(i=0,1,...,n). Gesucht ist die zugehörige Bézier-

Kurve.

1. Setze jj btbrr

=)(0 (j=0,1,...,n)

2. Berechne für r=1,...,n und j=0,...,n-r: 11

1 )()1()( −+

− +−= rj

rj

rj bttbttb

rrr

Dann ist )(0 tb nr

die Bézier-Kurve zu den gegebenen Kontrollpunkten.

BEISPIEL 3: MATLAB-PROGRAMM „BÉZIER-KURVEN“ MIT GRA FI-

SCHER OBERFLÄCHE

Wir schreiben ein MATLAB-Programm mit grafischer Oberfläche, welches dem Benut-

zer ermöglicht, die Kontrollpunkte per Mausklick einzugeben. Die Kontrollpunkte wer-

den durch einen Polygonzug („Kontrollpolygon“) verbunden, und es wird die zugehöri-

ge Bézier-Kurve gezeichnet. Klickt man auf einen schon vorhandenen Kontrollpunkt, so

kann man diesen an eine andere Stelle ziehen.

Lösung:


Algorithmus von de Casteljau: Datei bezier.m

function b=bezier(t,controlPoints)

%

% Berechnet die Bezier-Kurve zu den Kontrollpunkten p1,...,p(n+1)

% an der Stelle t

%

% Eingabeparameter: t Zeitpunkt

% controlPoints Matrix mit n+1 Kontrollpunkten

% Die Kontrollpunkte sind die Spalten der

% Matrix

%

% Rückgabewert: b Bezierkurve an der Stelle t

%

n = size(controlPoints); % 1. Index: 2 (ebene Kurve) oder 3 (räumliche Kurve)

% 2. Index: Anzahl Kontrollpunkte

numPoints = n(2)-1;

br = controlPoints;

bb = br; % bj0

for r=1:numPoints

br = zeros(n(1),numPoints-r+1);

for j=1:numPoints-r+1

br(:,j) = (1-t)*bb(:,j)+t*bb(:,j+1);

end

bb = br;

end

b=bb;

Erzeugen Sie ein leeres GUI. Ordnen Sie Achsen auf dem Fenster an, und stellen Sie XLimMode und YLimMode

auf manual ein. Wählen Sie als Einheit (Units) normalized. Ergänzen Sie Callbacks für das gesamte Fenster (View-

>Figure Callbacks) für die Ereignisse WindowButtonUp (Maustaste losgelassen) und WindowButtonMotion

(Mausbewegung).

function varargout = bezierDemobezierDemobezierDemobezierDemo(varargin)

% BEZIERDEMO M-file for bezierDemo.fig

% BEZIERDEMO, by itself, creates a new BEZIERDEMO or raises the existing

% singleton*.


%

% H = BEZIERDEMO returns the handle to a new BEZIERDEMO or the handle to

% the existing singleton*.

%

% BEZIERDEMO('CALLBACK',hObject,eventData,handles,...) calls the local

% function named CALLBACK in BEZIERDEMO.M with the given input arguments.

%

% BEZIERDEMO('Property','Value',...) creates a new BEZIERDEMO or

% raises the

% existing singleton*. Starting from the left, property value pairs are

% applied to the GUI before bezierDemo_OpeningFunction gets called. An

% unrecognized property name or invalid value makes property application

% stop. All inputs are passed to bezierDemo_OpeningFcn via varargin.

%

% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one

% instance to run (singleton)".

%

% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help bezierDemo

% Last Modified by GUIDE v2.5 06-Nov-2006 17:01:40

% Begin initialization code - DO NOT EDIT

gui_Singleton = 1;

gui_State = struct('gui_Name', mfilename, ...

'gui_Singleton', gui_Singleton, ...

'gui_OpeningFcn', @bezierDemo_OpeningFcn, ...

'gui_OutputFcn', @bezierDemo_OutputFcn, ...

'gui_LayoutFcn', [] , ...

'gui_Callback', []);

if nargin & isstr(varargin{1})

gui_State.gui_Callback = str2func(varargin{1});

end

if nargout

[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});

else

gui_mainfcn(gui_State, varargin{:});


end

% End initialization code - DO NOT EDIT

% --- Executes just before bezierDemo is made visible.

function bezierDemo_OpeningFcnbezierDemo_OpeningFcnbezierDemo_OpeningFcnbezierDemo_OpeningFcn(hObject, eventdata, handles, varargin)

% This function has no output args, see OutputFcn.




% varargin command line arguments to bezierDemo (see VARARGIN)

% Choose default command line output for bezierDemo

handles.output = hObject;

global points;

global numPoints;

numPoints = 0;

points= zeros(2,100);

global dragIndex;

dragIndex = -1;

% Update handles structure

guidata(hObject, handles);

% UIWAIT makes bezierDemo wait for user response (see UIRESUME)

% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.

function varargout = bezierDemo_OutputFcnbezierDemo_OutputFcnbezierDemo_OutputFcnbezierDemo_OutputFcn(hObject, eventdata, handles)

% varargout cell array for returning output args (see VARARGOUT);




% Get default command line output from handles structure

varargout{1} = handles.output;


% --- Executes on mouse press over figure background, over a disabled or

% --- inactive control, or over an axes background.

function figure1_WindowButtonDownFcnfigure1_WindowButtonDownFcnfigure1_WindowButtonDownFcnfigure1_WindowButtonDownFcn(hObject, eventdata, handles)

% hObject handle to figure1 (see GCBO)



global points;

global numPoints;

global dragIndex;

mouse = get(gca,'currentPoint');

p = [mouse(1,1);mouse(1,2)];

dragIndex = check(p, points,numPoints);

fprintf('dragIndex: %d',dragIndex);

if dragIndex == -1

numPoints = numPoints+1;

points(1,numPoints) = mouse(1,1);

points(2,numPoints) = mouse(1,2);

end

hold on;

c = get(gca,'children');


x = points(1,1:numPoints);

y = points(2,1:numPoints);

plot(x,y,'ro');

plot(x,y,'b');

controlPoints=points(:,1:numPoints);

t = 0:0.01:1;

y = zeros(2,101);

for i=1:101

y(:,i) = bezier(t(i),controlPoints);

end

plot(y(1,:),y(2,:),'g');

% --- Executes on mouse motion over figure - except title and menu.


function figure1_WindowButtonMotfigure1_WindowButtonMotfigure1_WindowButtonMotfigure1_WindowButtonMotionFcnionFcnionFcnionFcn(hObject, eventdata, handles)




global dragIndex;

global points;

global numPoints;

if dragIndex >= 1

c = get(gca,'children');


mouse = get(gca,'currentPoint');

p = [mouse(1,1);mouse(1,2)];

% fprintf('\nDrag: %f %f',p(1),p(2));

points(:,dragIndex)=p;

x = points(1,1:numPoints);

y = points(2,1:numPoints);

plot(x,y,'ro');

plot(x,y,'b');

controlPoints=points(:,1:numPoints);

t = 0:0.01:1;

y = zeros(2,101);

for i=1:101

y(:,i) = bezier(t(i),controlPoints);

end

plot(y(1,:),y(2,:),'g');

end

% --- Executes on mouse press over figure background, over a disabled or

% --- inactive control, or over an axes background.

function figure1_WindowButtonUpFcnfigure1_WindowButtonUpFcnfigure1_WindowButtonUpFcnfigure1_WindowButtonUpFcn(hObject, eventdata, handles)




global dragIndex;

dragIndex = -1;


Funktion zum Test, ob ein schon vorhandener Punkt geklickt wird: Datei check.m

function index = check= check= check= check(p,xySet, maxIndex)

%

% Prüft, ob der Punkt (x,y) in der

% Menge xySet enthalten ist

%

% Eingabeparameter: p: Koordinaten des zu prüfenden Punktes

% (Spaltenvektor)

% xySet: Matrix der Punkte. Ein Punkt ist als

% Spaltenvektor abgelegt

% maxIndex: Nur die Spalten 1...maxIndex von xySet

% werden durchsucht

%

% Rückgabewert: index des Punktes in xySet, falls Punkt gefunden,

% -1 sonst

n = size(xySet);

res = -1;

for i=1:maxIndex

dist = norm(p-xySet(:,i));

if dist < 0.05

res = i;

% fprintf('\nTreffer: %d',res);

break;

end

end

index = res;

Zieht man einen Punkt, so flackert das Fenster. Das Flackern wird abgestellt, wenn bei den figurefigurefigurefigure Eigenschaften

DoubleBufferingDoubleBufferingDoubleBufferingDoubleBuffering auf OnOnOnOn gesetzt wird.


Wie man sieht, sind Bézier-Kurven sehr gut kontrollierbar. Deshalb werden sie sie in

der Computergrafik häufig zur Konstruktion von Freiformkurven verwendet. Auch das

in WINDOWS enthaltene Zeichenprogramm PAINTBRUSH und das Zeichenpro-

gramm von OpenOffice unterstützen Bézier-Kurven.

Man kann zum Zweck der Interpolation auch mehrere Bézier-Kurven zusammensetzen.


3.6 ZWEIDIMENSIONALE I NTERPOLATION UND APPROXIMATION

3.6.1 ZWEIDIMENSIONALE LINEARE I NTERPOLATION

Die lineare Interpolation lässt sich leicht auf den zweidimensionalen Fall übertragen.

11.2

1.41.6

1.82

2

2.5

30

0.5

1

1.5

2

2.5

3

y x

(x0,y0,f(x0,y0))

(x1,y0,f(x1,y0)) (x1,y1,f(x1,y1))

(x,y,g(x,y))

Abbildung 9: Lineare Interpolation in zwei Variablen

Ist die Funktion f an drei Punkten ),(),,(),,( 110100 yxyxyx bekannt, so ergibt sich

der interpolierte Wert für [ ] [ ]1010 ,,),( yyxxyx ×∈ gemäß

01

01110

01

0001000

),(),()(

),(),()(),(),(

yy

yxfyxfyy

xx

yxfyxfxxyxfyxg

−−

−+−−

−+=

Dies ist die 3-Punkte-Form einer Ebene. Der Wert der linearen Interpolation am vierten

Punkt ist festgelegt. Daher ist diese Methode nicht geeignet zur Interpolation einer

Punktmenge, die auf einem Gitter gegebene ist.


3.6.2 BILINEARE I NTERPOLATION

Sind die Werte von f an vier Punkten ),(),,(),,(),,( 11100100 yxyxyxyx bekannt,

so wird bei der bilinearen Interpolation wie folgt vorgegangen, um einen Interpolati-

onswert an einer Stelle [ ] [ ]1010 ,,),( yyxxyx ×∈ zu berechnen:

1. Lineare Interpolation in x:

( ) ( )01

00010001

,),(),()(

xx

yxfyxfxxyxfxz

−−

−+=

( ) ( )01

10110102

,),(),()(

xx

yxfyxfxxyxfxz

−−

−+=

2. Lineare Interpolation in y:

( )01

1201

)()()(),(

yy

xzxzyyxzyxb

−−

−+=

Abbildung 10: Bilineare Interpolation


Insgesamt erhält man:

( ) ( )

( )( ) ( ) ( ) ( )

( )( ) ( )( ) ( )( ) ( )( )( )( )

( )( ) ( )( ) ( )( ) ( )( ) ( )( )( )( )

( )( ) ( )( ) ( )( ) ( )( )[ ]( )( )

( )( ) ( )( )[ ]( )( )

( )( ) ( )( )[ ]( )( )

( )( )[ ]( )( )

( )( )[ ] ( )( )[ ]( )( )

( )( )[ ] ( )( )[ ]( )( )0101

00111001

0101

10101100

0101

0011

0101

0001001

0101

0001010

0101

00001010010100

0101

00010000001101100

0101

100100001010010100

01

01

0001000

01

1011010

0

01

0001000

),(),(

),(),(

),(

),(),(

),(

,),(),(,),(

),(,),(),(

,),(),(

,),(),(

,),(),(),(

yyxx

yyxxyxfyyxxyxf

yyxx

xxyyyxfyyxxyxf

yyxx

yyxxyxf

yyxx

yyxxyyxxyxf

yyxx

yyxxxxyyyxf

yyxx

yyxxyyxxyyxxyyxxyxf

yyxx

yxfyxfyyxxyxfyyxxyxfyxfyyxx

yyxx

yxfxxyyyxfyxfyyxxyyxxyxf

yy

xx

yxfyxfxxyxf

xx

yxfyxfxxyxf

yy

xx

yxfyxfxxyxfyxb

−−−−+−−

+−−

−−+−−=

=−−

−−

+−−

−−−−−+

−−−−−−−

+−−

−−+−−−−−−−−=

−−−−−−−−−−−−

+−−

−−+−−−+−−=

−−−

−−−−−

−+−

+−−

−+=

3.6.3 ZWEIDIMENSIONALE BIKUBISCHE SPLINES

Problemstellung: Gesucht ist eine möglichst glatte Fläche, die durch eine vorgegebene

Menge von Punkten ),,( ijji uyx geht, wobei die ),( ji yx ein Rechteckgitter bilden:

dyyybxxxa mn =<<<==<<<= ...c ,... 1010


Abbildung 11: Rechteckgitter

Die bikubische Splinefunktion S = S(x,y) wird durch folgende Eigenschaften definiert

[ENGELN-MÜLLGES, REUTTER(1990)]:

1. ijji uyxS =),( (i=0,...,n, j=0,...,m)

2. S ist stetig differenzierbar

3. yx

S

∂∂∂ 2

ist stetig

4. In jedem Teilrechteck { }11,:),( ++ ≤≤≤≤= jjiiij yyyxxxyxR ist S

durch ein bikubisches Polynom ijp gegeben:

∑∑= =

−−=3

0

3

0

)()(),(k s

sj

kiijksij yyxxayxp

5. Zusätzlich erfüllt S noch bestimmte Randbedingungen. Häufig werden be-

stimmte partielle Ableitungen von S vorgegeben:


11

2

10

01

),(

),(

),(

ijijji

ijijji

ijijji

aryxyx

S

apyxx

S

aqyxy

S

==∂∂

∂

==∂∂

==∂∂

Durch die obigen Bedingungen ist die bikubische Splinefunktion eindeutig bestimmt.

Die MATLAB-Funktion interp2 führt eine bikubische Interpolation durch:

Grobes Gitter >> [x,y]=meshgrid(-3:1:3);

Beispiel-Funktion auf

grobem Gitter >> z=peaks(x,y);

Feines Gitter >> [xFine,yFine]=meshgrid(-3:.25:3);

Lineare Interpolation

auf feinem Gitter

>> zLin = interp2(x,y,z,xFine,yFine,'linear');

>> mesh(x,y,z);

>> hold on;

>> mesh(xFine,yFine,zLin+10);


Abbildung 12: Lineare Interpolation

Bikubische Interpola-

tion auf feinem Gitter

>> zBicubic = interp2(x,y,z,xFine,yFine,'spline');

>> mesh(x,y,z);

>> hold on;

>> mesh(xFine,yFine,zBicubic+10);


Abbildung 13: Bikubische Splineinterpolation

3.6.4 ANWENDUNG : VERGRÖßERN DIGITALER BILDER

Bild laden

orig = imread('test.jpg');

orig = double(orig); % in double-Werte umwandeln

size(orig)

ans =

670 1000 3

% Erster Index : Anzahl Zeilen

% Zweiter Index: Anzahl Spalten

% Dritter Index: 1: Rot 2: Grün 3: Blau

image(orig/255); % Bei double-Bild: Farben zwischen 0 und 1


Abbildung 14: Original-Bild

Gitter, auf

dem das Ori-

ginal defi-

niert ist:

breite = size(orig,2)

breite =

1000

hoehe = size(orig,1)

hoehe =

670

[xOrig,yOrig]=meshgrid(1:breite,1:hoehe);

Bereich, der

vergrößert

werden soll:

[xZoom,yZoom]=meshgrid(350:.2:450,350:.2:450);

Zoom mit

Interpolati-

onsmethode

‚nearest’

% Jede Farbe einzeln

z1(:,:,1)=interp2(xOrig,yOrig,orig(:,:,1),xZoom,yZoom,'nearest');




(stückweise

konstant)

image(z1/255);

>>

Abbildung 15: Vergrößerung mit Interpolationsmethode 'nearest'

Zoom mit

Interpolati-

onsmethode

‚spline’

z2(:,:,1)=interp2(xOrig,yOrig,orig(:,:,1),xZoom,yZoom,'spline');



% Hier ist nicht garantiert, dass alle Werte zwischen 0 und 255

% liegen

totalMin = min(min(min(z2)));

z2=z2-totalMin;

totalMax = max(max(max(z2))) ;

image(z2/totalMax)


Abbildung 16: Vergrößerung mit Interpolationsmethode 'Spline'

Professionelle Bildbearbeitungsprogramme unterstützen bikubische Splines. In man-

chen Fällen findet man weitere Spezialvarianten wie sog. S-Splines.


AUFGABE SB4-4: BIKUBISCHE/BILINEARE INTERPOLATION

Erzeugen Sie eine Anwendung mit grafischer Oberfläche, welche die bikubische und die

bilineare Interpolation demonstriert.

Mit zwei Radio Buttons soll zwischen bilinearer und bikubischer Interpolation umge-

schaltet werden können.

Als Ausgangspunkt soll die Ebene durch die 16 Punkte ( ) )0,1,1(0,, −−= jiyx ji

(i=1,...,4, j=1,...,4) dienen.

Ihr Programm soll mit 16 Radio Buttons die Auswahl des Punktes, dessen z-Koordinate

verändert werden soll, ermöglichen. Die z-Koordinate soll mit einem Schieberegler

zwischen -5 und +5 verändert werden können. Wird der Schieberegler betätigt, so soll

sofort die interpolierte Fläche auf dem Gitter ( )25.0*,25.0* ji (i,j=0,...,12) geplottet

werden.

Anleitung: Leider unterstützt MATLAB 6.5 noch nicht das Konzept der Button Group,

welches es erlaubt, mehrere Radio Buttons logisch zu einer Gruppe zusammenzufassen,

so dass bei Drücken eines Buttons automatisch alle anderen Buttons der Gruppe in den

Zustand „AUS“ gehen. Dieses Verhalten muss bei MATLAB 6.5 noch „zu Fuß“ pro-

grammiert werden.

Wenn Sie also MATLAB 6.5 benutzen, gehen Sie wie folgt vor:

Legen Sie dazu in der OpeningFcn eine globale Variable an, welche die Handles der

Radio Buttons enthält:

global hRadioButtons;

hRadioButtons = zeros(1,16);

In der Callback-Funktion jedes Radio Button wird das zugehörige Handle beim ersten

Aufruf initialisiert:

function radiobutton1_Callback(hObject, eventdata, function radiobutton1_Callback(hObject, eventdata, function radiobutton1_Callback(hObject, eventdata, function radiobutton1_Callback(hObject, eventdata,

handles)handles)handles)handles)

��

Aufgabe

��

20 Punkte


% hObject handle to radiobutton1 (see GCBO)

% eventdata reserved - to be defined in a future

version of MATLAB

% handles structure with handles and user data

(see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of

radiobutton1

global hRadioButtons;

if hRadioButtons(1) == 0

hRadioButtons(1) = hObject;

end

Rufen Sie ferner in der Callback-Funktion eines Buttons folgende Funktion auf:

update(hRadioButtons,i); % i: Index des Radio Button

Dabei setzt die Funktion update (Datei update.m) alle Buttons außer einem auf „AUS“:

function update(hRadioButtons,i)

% Alle Buttons außer dem i-ten auf “AUS” setzen

for k=1:16

if hRadioButtons(k) == 0

continue;

end

set(hRadioButtons(k),'Value',0);

end

set(hRadioButtons(i),'Value',1);


4 Z USAMME NFASSUN G

In diesem Kapitel wurden verschiedene Verfahren der Interpolation vorgestellt. Bei der

Interpolation wird stets verlangt, dass die konstruierte Funktion durch alle vorgegebe-

nen Punkte geht. Lineare Interpolation ist die einfachste Variante. An den Stützstellen

ist die stückweise lineare Interpolation nicht differenzierbar. Zur Konstruktion glatter

Kurven greift man daher oft auf andere Verfahren zurück. Die Polynominterpolation

liefert zwar glatte Funktionen, neigt aber mit wachsender Anzahl Stützstellen zu starken

Oszillationen. Daher werden in der Praxis häufig kubische Splines eingesetzt. Diese

setzen sich aus Polynomen vom Grad 3 zusammen und kommen damit dem, was man

intuitiv unter einer möglichst glatten interpolierenden Kurve versteht, sehr nahe.

Das Konzept der kubischen Splines lässt sich auch auf Kurven in Ebene und Raum

sowie auf Flächen übertragen.

Bézier-Kurven sind ein Beispiel für die Konstruktion von Freihand-Kurven. Hier geht

es weniger um Interpolation als darum, die Form einer Kurven möglichst gut kontrollie-

ren zu können.


5 L I TE R AT UR

Schweizer, J.: Numerik für Bioinformatiker. 2002

http://homepages.uni-tuebingen.de/juergen.schweizer/nfb.html

(zuletzt besucht am 27.10.2006)

MathWorks: Building GUIs with MATLAB. 2006.

www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/buildgui.pdf


Deutsche Welle: Achterbahn. Projekt Zukunft.

http://www.deutsche-welle.de/dw/article/0,2144,1789282,00.html


Babovsky, H., Neundorf, W.: Numerische Approximation von Funktionen

(Mathematische Methoden des Computer Aided Geometric Design). TU Ilme-

nau, 2006.

http://www3.tu-ilmenau.de/

/site/phys/fileadmin/template/ifm/user/Neundorf/Dokumente/Lehre/num_app/n

um_app.pdf


Engeln-Müllges, G., Reutter, F.: Formelsammlung zur Numerischen Mathema-

tik mit C-Programmen. BI Wissenschaftsverlag, 2. Aufl., 1990.

Sormann, H.: Numerische Methoden in der Physik. 2006

http://itp.tugraz.at/LV/sormann/NumPhysik/Skriptum/kapitel3.pdf


Verfürth, R.: Einführung in die Numerische Mathematik. Vorlesungsskriptum Sommer-

semester 2005. Ruhr-Universität Bochum.

http://www.ruhr-uni-bochum.de/num1/skripten/Numerik_Einfuehrung.pdf


http://homepages.uni-tuebingen.de/juergen.schweizer/nfb.html

http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/buildgui.pdf

http://www.deutsche-welle.de/dw/article/0,2144,1789282,00.html

http://www3.tu-ilmenau.de/site/phys/fileadmin/template/ifm/user/Neundorf/Dokumente/Lehre/num_app/num_app.pdf



http://itp.tugraz.at/LV/sormann/NumPhysik/Skriptum/kapitel3.pdf

http://www.ruhr-uni-bochum.de/num1/skripten/Numerik_Einfuehrung.pdf

SB4 Interpolation

Documents

Transcript of SB4 Interpolation