%% Blatt 8 - Programmieraufgaben:
%% Splines nach Mansfield, Cox & de Boor (4 Punkte)
% Vorname Nachname Mat.-Nr.
% Wir wollen in diesen Programmieraufgaben mit den allgemeinen Basis Splines, den B-Splines auf nicht-äquidistanten Gittern, rechnen. Die Benutzung von Splines in dieser Darstellung ist in der Computergraphik und den CAx-Technologien weit verbreitet. Dazu betrachten wir das Knoten-Gitter
% \Delta = (t_1, t_2, \ldots, t_{n+1})  mit  t_1 \leq t_2 \leq \cdots \leq t_{n+1}.
% Die B-Splines vom Grad 0 sind dann durch
% B_{0, j}(x) = \begin{cases} 1 , & \text{falls } t_j \leq x < t_{j+1} , \\ 0 , & \text{sonst,} \end{cases}
% mit j=1, 2, \ldots, n gegeben. Für die Splines von höherem Grad erweitern wir das Knotengitter durch Repetition des Anfangs- und Endknotens, das heisst wir setzen
% \hat{t}_j := t_{\min\big(n+1, \max(1, j)\big)} für jedes j \in \mathbb{Z}.
% Gemäss der Cox-de Boor Rekursionsformel ergeben sich dann die B-Splines höheren Grades durch die Rekursionsformeln
% B_{m, 1}(x) = [1 - w_{m-1, 1}(x)] B_{m-1, 1}(x)  \\  B_{m, j}(x) = w_{m-1, j-1}(x) B_{m-1, j-1}(x) + [1 - w_{m-1, j}(x)] B_{m-1, j}(x)  \\  B_{m, n+m}(x) = w_{m-1, n+m-1}(x) B_{m-1, n+m-1}(x)
% mit m \geq 1 und j = 2, \ldots, n+m-1, dabei ist
% w_{m, j}(x) = \begin{cases} \frac{x - \hat{t}_j}{\hat{t}_{j+m} - \hat{t}_j} , & \text{falls } \hat{t}_j \neq \hat{t}_{j+m} , \\ 0 , & \text{sonst.} \end{cases}
% Da die Rekursion in den Rekursionsformeln eine Dreiecksstruktur aufweist, kann man die Rekursion auflösen, um auf folgenden Algorithmus fuer die Berechnung von B_{m,1}(x), \ldots, B_{m, n+m}(x) zu gelangen:
% Input: Gitter \mathbf{t} = [t_{j}]_{j=1}^{n+1} \in \mathbb{R}^{n+1}, Grad m \in \mathbb{N}, Auswertungsorte \mathbf{x} = [x_{i}]_{i=1}^{\ell} \in \mathbb{R}^{\ell}
% Output: Auswertungswerte \mathbf{y} = [B_{m, j}(x_i)]_{i,j=1}^{\ell, n+m} \in \mathbb{R}^{\ell \times (n+m)}
% initialisiere \mathbf{y} := \mathbf{0} \in \mathbb{R}^{\ell \times (n+m)}
% setze y_{i, j} := 1 für alle i=1, \ldots \ell und j = 1,\ldots,n mit t_j \leq x_i < t_{j+1}
% for k = 1, \ldots , m do
% \qquadfor j = n+k-1, n+k-2, \ldots, 1 do
% \qquad\qquadsetze tl := t_{\min\big(n+1, \max(1, j-k+1)\big)}
% \qquad\qquadsetze tr := t_{\min\big(n+1, \max(1, j+1)\big)}
% \qquad\qquadif tl \neq tr
% \qquad\qquad\qquadsetze \mathbf{w} = (tr - \mathbf{x}) / (tr - tl)
% \qquad\qquad\qquadsetze \mathbf{y}_{:, j+1} = \mathbf{y}_{:, j+1} + (1 - \mathbf{w}) * \mathbf{y}_{:, j}
% \qquad\qquad\qquadsetze \mathbf{y}_{:, j} = \mathbf{w} * \mathbf{y}_{:, j}
% \qquad\qquadelse
% \qquad\qquad\qquadsetze \mathbf{y}_{:, j+1} = \mathbf{y}_{:, j+1} + \mathbf{y}_{:, j}
% \qquad\qquad\qquadsetze \mathbf{y}_{:, j} = \mathbf{0}
% \qquad\qquadend if
% \qquadend for
% end for
% Der Operator * steht hierbei für die elementweise Multiplikation zweier Vektoren.
% Schreiben Sie hier am Ende dieses Live Scripts eine MATLAB-Funktion,
% ys = bSplines(ts, m, xs) ,
% welche den obigen Algorithmus umsetzt; dabei sind die Ein- und Ausgaben der Form
% \mathtt{ts} = \pmatrix{ t_1 \cr t_2 \cr \vdots \cr t_{n+1} } , \mathtt{xs} = \pmatrix{ x_1 \cr x_2 \cr \vdots \cr x_{\ell} }  und  \mathtt{ys} = \pmatrix{B_{m,1}(x_1) & \cdots & B_{m,n+m}(x_1) \cr B_{m,1}(x_2) & \cdots & B_{m,n+m}(x_2) \cr \vdots & \ddots & \vdots \cr B_{m,1}(x_\ell) & \cdots & B_{m,n+m}(x_\ell)}.
%
% Der nachfolgende Code zeichnet die B-Splines auf einem Gitter mit n=5 für die Grade m = 0, 1, 2, 3:

ts = ((0:5) + .6*rand(1, 6)).';
xs = linspace(min(ts)-0.2, max(ts)+0.2, 2^11+1).';
figure();
p = 0;
tic;
for m=0:1:3
    p = p+1;
    subplot(2, 2, p);
    ys = bSplines(ts, m, xs);
    plot(xs, ys);
end
toc;

% Ausgehend davon lassen sich nun parametrische Kurven zeichnen. Für m \in \mathbb{N} betrachtet man die Kontrollpunkte \mathbf{p}_1, \ldots, \mathbf{p}_{n+m} \in \mathbb{R}^2, die man den Kontrollpolygonzug nennt, und definiert den Spline vermittels
% s(t) := \sum_{j=1}^{n+m} \mathbf{p}_j B_{m, j}(t)
% für t \in [t_1, t_{n+1}]. In nachfolgendem Bild sind für j=1,\ldots,n+1 einfach t_j = j-1 gewählt. Wenn Sie die Check-Box bei der Variable extern aktivieren, wird ein externes Graphikfenster geöffnet, worin Sie mit Links-Click&Drag die Kontrollpunkte verschieben und mit Rechts-Klick in das Fenster durch die verschiedenen, möglichen m durch iterieren können.

extern = false;
m = 3;
xc = [.1; 1; 1.6; .2; 1.6; 1; .1];
yc = [.2; 2; .4; 0; -.4; -3; -.2];
splineFigure(m, xc, yc);
if extern
    set(gcf,'Visible','on');
end

%
% Funktionen:

function ys = bSplines(ts, m, xs)
    %% TODO
    
end

% Funktionem für interaktive Splines:

function splineFigure(m, xc, yc)
    fH = figure();
    set(fH, 'WindowButtonDownFcn', @MouseClick);
    set(fH, 'WindowButtonUpFcn', @MouseRelease);
    aH = axes();
    data.m = m;
    data.xc = xc;
    data.yc = yc;
    updateSplineData();
    [pH, tH] = initSplineDraw();

    function updateSplineData()
        data.l = size(data.xc, 1);
        data.n = data.l-data.m;
        data.tc = (0:data.n).';
        data.ts = linspace(0, (1-eps)*data.n, 2^11+1).';
        data.bs = bSplines(data.tc, data.m, data.ts);
        data.xs = data.bs * data.xc;
        data.ys = data.bs * data.yc;
    end
    function [pH, tH] = initSplineDraw()
        pH = plot(data.xc, data.yc, 'k:s', data.xs, data.ys, 'r-');
        axis('equal');
        tH = title(['Spline mit Grad m = ' num2str(data.m)]);
    end
    function updateSplineDraw()
        set(pH(1), 'XData', data.xc, 'YData', data.yc);
        set(pH(2), 'XData', data.xs, 'YData', data.ys);
        set(tH, 'String', ['Spline mit Grad m = ' num2str(data.m)]);
        drawnow
    end
    function moveSplineData(x, y)
        data.xc(data.i) = x;
        data.yc(data.i) = y;
    end
    function MouseClick(src, evnt)
        if strcmp(get(src,'SelectionType'),'normal')
            cp = get(aH, 'CurrentPoint');
            x = cp(1, 1);
            y = cp(1, 2);
            [~, data.i] = min((data.xc-x).^2 + (data.yc-y).^2);
            moveSplineData(x, y);
            updateSplineData();
            updateSplineDraw();
            set(src,'WindowButtonMotionFcn', @MouseHold);
        elseif strcmp(get(src,'SelectionType'),'alt')
            if data.m < data.l-1
                data.m = 1 + data.m;
            else
                data.m = 1;
            end
            updateSplineData();
            updateSplineDraw();
        end
    end
    function MouseHold(src, evnt)
        cp = get(aH, 'CurrentPoint');
        moveSplineData(cp(1, 1), cp(1, 2));
        updateSplineData();
        updateSplineDraw();
    end
    function MouseRelease(src, evnt)
        set(src,'WindowButtonMotionFcn','');
    end
end