function doIt = simtrans2()
    clear all

    global w h wn hn sn T P imgQn imgDn Gdx Gdy
    w = 10;
    h = 10;

    wn = w - 1;
    hn = h - 1;
    sn = max(wn,hn);

    %Parameter init
    t = 0;
    T = [1, 0, t; 0, 1, t];
    P = [t, t];
    %Ziel Transformation
    targetT = [1, 0, -3; 0, 1, -3];

    %Bilder erstellen
    imgQ(1:w*h) = 128;
    imgQ = reshape(imgQ,h,w);
    imgQ(2:4,2:4) = 255;
    imgD = transImg(imgQ, targetT);

    %Bilder normieren
    imgQn = imgQ / 255.;
    imgDn = imgD / 255.;

    %Gradienten erstelle
    [Gdx, Gdy] = gradient(imgDn);

    doIt=@doIteration;
    fprintf('simtrans initialized\n')


%normierte Koordinaten in Matrixkoordinaten umrechnen
function [xp, yp] = deNormCoords(xn,yn)
	global wn hn sn
	xp = (xn * sn + wn) / 2 + 1; %Matrix Index beginnt bei 1
	yp = (yn * sn + hn) / 2 + 1; %Matrix Index beginnt bei 1


%Koordinaten um T verschieben
function [xt, yt] = transCoords(xn, yn, T)
	V = [xn, yn, 1];
	R = T * V';
	xt = R(1,1);
	yt = R(2,1);


%Punkt der Matrix M an der Stelle x,y bilinear interpolieren
function v = biInt(M, x, y)
	global w h

	if ((x < 1) || (y < 1) || (x > w) || (y > h)) 
		v = 0;
	else
		v = interp2(M, x, y);
    end


%Ein Bild um T transformieren
function [It] = transImg(img, T)
	global wn hn
	
	stepX = 2 / wn;
	stepY = 2 / hn;
	for xn=-1:stepX:1
		for yn=-1:stepY:1
			[xp, yp]   = deNormCoords(xn,yn);				%Nicht transformierte Koordinaten denormalisieren										
			[xpt, ypt] = transCoords(xp, yp, T);				%Koordinaten mit T transformieren
			It(round(xp), round(yp)) = biInt(img, xpt, ypt);	
        end
    end		


%Jacobi Matrix an der Stelle x,y erstellen
function JC = jacobi(x,y)
	JC = [1,0;0,1];


%Ein Iterationsschritt
function deltaP = doIteration()
	global wn hn Gdx Gdy imgQn imgDn P T
	%imgQn ist das Quellbild normiert auf Werte zwischen 0 und 1
	%imgDn ist das Zielbild normiert auf Werte zwischen 0 und 1
	%P ist der Parameter Vektor, welcher bei der Translation (tx,ty) ist
	%T ist die Transformationsmatrix = [1,0,tx;0,1,ty]
	%Gdx Gdy entpricht gradient(imgDn);

	A = zeros(2,2);
	b = zeros(2,1);
	stepX = 2 / wn;							%2 / (Bildbreite-1)
	stepY = 2 / hn;							%2 / (Bildh??he-1)

	for xn=-1:stepX:1							%ber die normierten Pixelkoordinaten laufen
		for yn=-1:stepY:1
			[xnt, ynt] = transCoords(xn, yn, T);			%Transformationsmatrix auf normierte Koordinaten anwenden
			[xp, yp]   = deNormCoords(xn,yn);			%Nicht transformierte Koordinaten denormalisieren							
			[xpt, ypt] = deNormCoords(xnt, ynt);			%transformierte Koordinaten denormalisieren			
			JC = jacobi(xn, yn);					%Jacobi Matrix erzeugen
			GD = [biInt(Gdx, xpt, ypt), biInt(Gdy, xpt, ypt)];	%Gradienten Vektor vom ZielBild erzeugen, biInt interpoliert den Wert bilinear, 0 wenn x,y au??erhalb
			
			A = A + JC' * (GD' * GD) * JC;
			b = b + JC' * (biInt(imgDn, xpt, ypt) - biInt(imgQn, xp, yp)) * GD'; %biInt interpoliert den Wert bilinear, 0 wenn x,y au??erhalb

        end
    end
	A
	b = -b
	deltaP = A \ b;
	P = P + deltaP';

	T = [1, 0, P(1,1); 0, 1, P(1,2)]