Verfasst am: 21.08.2012, 15:04
Titel: Vektorisieren von for-Schleifen
Hallo,
ich habe ein paar sehr zeitaufwendige for-Schleifen in meinem Programm laufen, und möchte das Programm gerne beschleunigen. Dabei habe ich an Vektorisieren gedacht, allerdings habe ich das noch nicht hinbekommen.
Q = [R(yN1,xN1),R(yN1,xN2);R(yN2,xN1),R(yN2,xN2)]; %Matrix mit den Werten der Stuetzpunkte
p = [(xN2-xN(i,j))/(xN2-xN1);(xN(i,j)-xN1)/(xN2-xN1)]; %Vektor mit Gewichtungsfunktionen fuer x-Richt.
o = [(yN2-yN(i,j))/(yN2-yN1),(yN(i,j)-yN1)/(yN2-yN1)]; %Vektor mit Gewichtungsfunktionen fuer y-Richt.
H(round(yO(j)),round(xO(i))) = o*(Q*p); %Bilineare Interpolation in Matrixschreibweise
end end end end
und:
for k=1:size(CoordBrake,1) if CoordBrake(k,1) < sizeX && CoordBrake(k,2) < sizeY && CoordBrake(k,1) > 0 && CoordBrake(k,2) > 0
Matrix(CoordBrake(k,2),CoordBrake(k,1)) = Replace(CoordBrake(k,2),CoordBrake(k,1));
end end
Wie gesagt, was ich bisher probiert habe hat nicht funktioniert, wäre für Hilfe sehr dankbar. Vielleicht weiß jemand auch eine andere Möglichkeit diesen Code zu beschleunigen.
ich habs leider immer noch nicht hinbekommen. Zwar konnte ich durch ändern der Reihenfolge der for-Schleifen und vordefinieren der sinus- und cosinus-Werte diese eine geschachtelte Schleife etwas beschleunigen, bei der Vektorisierung bin ich jedoch noch nicht weiter gekommen. Ich frage mal anders:
for j=1:length(y) %Hier Drehrichtung im Uhrzeigersinn!
for i=1:length(x)
xN(i,j) = round(x(i)*co+y(j)*si+origin(1)); %minus sinus in dieser Zeile: Gegenuhrzeigersinn
yN(i,j) = round(y(j)*co-x(i)*si+origin(2)); %minus sinus in dieser Zeile: Uhrzeigersinn
end end
Ich habe eine Matrix, die ich um einen Winkel theta rotieren möchte. Jetzt brauche die Koordinaten der neuen (rotierten) Matrix.
xN und yN sind ähnlich wie beim Befehl meshgrid zwei Matrizen die einem alten Koordinatenpaar ein neues (rotiertes) Paar zuordnen.Die oben angeführte Schleife liefert diese Koordinaten auch korrekt allerdings würde ich das gerne beschleunigen.
Ist hier eine Vektorisierung in dieser Form überhaupt möglich? Könnte man diese Rotation effizienter gestalten?
D.h. eine Matrix ist die Abbildungsmatrix, die man mit den Daten einmal nur initialisieren muss und dann kann man mit Hilfe der entsprechenden Zugriffsoperatoren auf die Spalten / Zeilen der Datenmatrix zugreifen und die Daten der Matrix abbilden
ich rotiere auch nicht die Matrix, da hab ich mich schlecht ausgedrückt, sondern das Koordinatensystem. Ich habe eine dreidim. Matrix mit 240x320x3000 Einträgen. Jede dieser 240x320 Matrizen entspricht einem Bild, welches mit der Zeit rotiert. Die Winkel sind mir bekannt, auch gelingt es mir mit dem obigen Skript die Koordinaten im neuen Bild zu ermitteln (d.h, das Koordinatensystem so zu drehen, dass ich aus dem ehemals rotierendem Bild ein stehendes erzeugen kann): Ich rotiere das Koordinatensystem greife die Werte an den neuen Positionen ab und schreibe sie auf die Positionen im nicht-rotierten Koordinatensystem.
Das Skript erfüllt also seinen Zweck, ich hätte nur gerne gewusst, ob dies auch schneller, effizienter geht.
Jede Spalte (Bild) Deiner Matrix ist ein Vektor und wenn Du die Winkel kennst, dann kannst Du es auch rotieren. Du wirst ja das ganze Bild rotieren und nicht nur einen Teilbereich, d.h. jedes Pixel (Element des Vektors) wird um den gleichen Winkel verändert.
D.h. Du musst die Rotationsmatrix passend bestimmen und dann auf die Datenmatrix anwenden (ggf, wenn sich die Matrix von Vektor zu Vektor bzw. Bild zu Bild ändert, brauchst Du eine Schleife)
Matematisch
Welches ist nun der Code, den Du beschleunigen möchtest?
In dem Code, den Du zeigst, werden die Ausgaben nicht pre-alloziert. Dies wäre ein dramatischer Nachteil. Beispiel:
Ich habe also einfach nur den Index aus der "for i"-Schleife in die Berechnung einbezogen. Dann kann man natürlich auch statt "x(1:length(x))" einach "x" schreiben. Eventuell muss man jedoch x noch transponieren.
Und die zweite Schleife läßt sich auch noch vermeiden:
Dann ist die Pre-allocation nicht hilfreich. Eventuell müssen auch hier wieder x und y anders transponiert werden. Aber das kannst Du am besten selbst testen.
Es ist immer hilfreich, genau den Code zu posten, der optimiert werden soll und passende Test-Daten mitzuliefern, damit wir die direkt laufen lassen können.
Q = [R(yN1,xN1),R(yN1,xN2);R(yN2,xN1),R(yN2,xN2)]; %Matrix mit den Werten der Stuetzpunkte
p = [(xN2-xN(i,j))/(xN2-xN1);(xN(i,j)-xN1)/(xN2-xN1)]; %Vektor mit Gewichtungsfunktionen fuer x-Richt.
o = [(yN2-yN(i,j))/(yN2-yN1),(yN(i,j)-yN1)/(yN2-yN1)]; %Vektor mit Gewichtungsfunktionen fuer y-Richt.
H(round(yO(j)),round(xO(i))) = o*(Q*p); %Bilineare Interpolation in Matrixschreibweise
end end end end
Dabei werden die Daten von R mit Hilfe der rotierten Koordinaten in xN und yN in eine neue Matrix H geschrieben, dabei wird noch eine bilineare Interpolation durchgeführt. Ich weiß Matlab bietet eine eigene Funktion dafür, konnte diese aber noch nicht erfolgreich einbauen, also hab ich es so programmiert.
Wäre super wenn du auch hierfür ein paar Tips zum Beschleunigen hättest. Vielen Dank schon mal und freundliche Grüße,
Die Definition von "origin" muss vor der Benutzung erfolgen.
Mit den Testdaten wird die eigentliche Berechnung nicht erreicht und H nicht definiert. Kannst Du etwas passenderes angeben?
Ein paar Dinge lassen sich trotzdem finden:
Code:
match = round(yO) < sizeY && round(yO) > 2; % Nur einmal testen! for i=1:lx
ifround(xO(i)) < sizeX && round(xO(i)) > 2% Ausserhalb der Schleife! for j=1:ly
if match(j) % yN und xN sind bereits gerundet! if yN(i,j) < sizeY && xN(i,j) < sizeX && yN(i,j) > 2 && xN(i,j) > 2 ...
vielen Dank für deine Hilfe, ich lerne dabei gerade wieder eine Menge dazu.
Leider konnte ich keine Dateien anhängen (funktionierte bei mir aus irgendeinem Grund nicht), deshalb habe ich eine Testmatrix im richtigen Format mit Zahlen zwischen 200 und 800 in dieser Schleife erzeugt:
Code:
M = zeros(240,320);
for i=1:320 for j=1:240
r = 200+(800-200).*rand(1,1);
M(j,i) = r;
end end
Q = [R(yN1,xN1),R(yN1,xN2);R(yN2,xN1),R(yN2,xN2)]; %Matrix mit den Werten der Stuetzpunkte
p = [(xN2-xN(i,j))/(xN2-xN1);(xN(i,j)-xN1)/(xN2-xN1)]; %Vektor mit Gewichtungsfunktionen fuer x-Richt.
o = [(yN2-yN(i,j))/(yN2-yN1),(yN(i,j)-yN1)/(yN2-yN1)]; %Vektor mit Gewichtungsfunktionen fuer y-Richt.
H1(round(yO(j)),round(xO(i))) = o*(Q*p); %Bilineare Interpolation in Matrixschreibweise
end end end end
Habe deinen Code schon ausprobiert, er läuft mit der oben angeführten Testmatrix auch etwas schneller, denke dass das meiste Potetial aber in der Berechnung im Schleifenkörper steckt.
habe versucht diese bipolare Interpolation mit der Matlab-Funktion interp2 umzusetzen, was mir aber nicht gelungenist:
ZI = interp2(X,Y,Z,XI,YI)
X und Y enthalten die Koordinaten der Stützstellen, was bei mir xN und yN entsprechen sollte, Z die oben angegebene Matrix. Ich scheitere leider schon daran xN und yN in ein passendes Format zu wandeln, sodass ich interp2 darauf anwenden kann, dann müssen diese Koordinaten in xN und yN ja noch mit ceil bzw. floor bearbeitet werden (wie oben im Code).
Dann habe ich da noch eine weitere Frage zur Vektorisierung:
Code:
rIn = 10;
zeta = linspace(0,pi/2,100);
origin = [154120];
Matrix = abs(randn(500,500)*100);
xR = round((rIn+10)*cos(zeta)+origin(1));
yR = round((rIn+10)*sin(zeta)+origin(2));
for j=1:size(zeta,2)
celsius2(j) = Matrix(yR(j),xR(j));
end
Was geschieht eigentlich, wenn xN(i,j) ein Integer ist? Dann ist "floor(xN(i,j)) == ceil(xN(i,j))". Ist das so gewollt?
Code:
origin = [154, 120];
R = M;
[sizeY, sizeX] = size(R); % TRANSPOSE costs time
x = -200:0.5:200;
y = -200:0.5:200;
xO = x + origin(1);
yO = y + origin(2);
lx = length(x);
ly = length(y);
theta = 0.1 * pi;
co = cos(theta);
si = sin(theta);
% "xN und yN müssen hier nicht gerundet werden" % Das sollten sie aber! Ansonsten muss dies wiederholt in der Schleife % gemacht werden!
xN = bsxfun(@plus, x'*co, y * si) + origin(1);
yN = bsxfun(@minus, y*co, x' * si) + origin(2);
Q = [R(yN1,xN1),R(yN1,xN2);R(yN2,xN1),R(yN2,xN2)]; %Matrix mit den Werten der Stuetzpunkte
p = [(xN2-xN(i,j))/(xN2-xN1);(xN(i,j)-xN1)/(xN2-xN1)]; %Vektor mit Gewichtungsfunktionen fuer x-Richt.
o = [(yN2-yN(i,j))/(yN2-yN1),(yN(i,j)-yN1)/(yN2-yN1)]; %Vektor mit Gewichtungsfunktionen fuer y-Richt.
H1(r_yO(j), r_xO(i)) = o*(Q*p); %Bilineare Interpolation in Matrixschreibweise
end end end end
Dann würde ich noch versuchen "o*Q*p" elementweise auszudrücken. Das Erstellen der kleinen temporären Vektoren benötigt einen unangemessen großen Overhead und der Aufruf einer BLAS-Funktion für so kleine Vektoren und Matrizen ist ebenfalls ineffizient.
Bitte ausgiebig testen, da ich das nur im Forum-Editor zusammengestückelt habe.
Meine Strategie ist also zunächstmal das De-Vektorisieren, also das genau Gegenteil von dem Gewünschten. Ich gehe aber davon aus, dass dies bereits schneller ist, da Matlab Schleifen sehr effizient bearbeiten kann. Wieviel bringt dies?
nocheinmal recht herzlichen Dank an dieser Stelle für deine Hilfe! Ich habe deinen Code jetzt ausprobiert, musste dabei ein paar Sachen ändern: Erstens, für unser Problem eigentlich unwesentlich, habe ich zwei Einträge in der Matrix Q (bzw. dann auch in der ausgeschriebenen Form) vertauscht. Zweitens: Danke für den Hinweis, wenn xN(i,j) (oder yN(i,j)) ein Integer ist sind ceil und floor gleich, d.h. aber das der Punkt für den interpoliert werden soll mit einem bekannten Datenpunkt zusammenfällt (es muss also gar nicht interpoliert werden, man kann den Punkt so nehmen wie er ist). Das hat in der ursprünglichen Version zu NaN-Einträgen geführt! Jetzt wird das einfach durch eine if-Abfrage verhindert.
Ich poste hier noch einmal verschiedene Versionen des Codes:
Code:
M = zeros(240,320);
for i=1:320 for j=1:240
r = 200+(800-200).*rand(1,1);
M(j,i) = r;
end end
origin = [154120];
[sizeX sizeY] = size(M');
x = -200:0.5:200;
y = -200:0.5:200;
xO = x+origin(1);
yO = y+origin(2);
lx = length(x);
ly = length(y);
theta = 0.1*pi;
co = cos(theta);
si = sin(theta);
xN = bsxfun(@plus,x'*co,y*si)+origin(1);
yN = bsxfun(@minus,y*co,x'*si)+origin(2);
tic
H1 = zeros(240,320);
for i=1:lx
for j=1:ly
ifround(yO(j)) < sizeY && round(xO(i)) < sizeX && round(yO(j)) > 1 && round(xO(i)) > 1 if yN(i,j) < sizeY && xN(i,j) < sizeX && yN(i,j) > 1 && xN(i,j) > 1
%______________________________________________________
%Bilineare Interpolation:
xN1 = floor(xN(i,j)); %x-Koord. des1. Stuetzpunktes
xN2 = ceil(xN(i,j)); %x-Koord. des2. Stuetzpunktes
yN1 = floor(yN(i,j)); %y-Koord. des1. Stuetzpunktes
yN2 = ceil(yN(i,j)); %y-Koord. des2. Stuetzpunktes if xN1 == xN2 || yN1 == yN2
H1(round(yO(j)),round(xO(i))) = M(yN1,xN1);
else
Q = [M(yN1,xN1),M(yN2,xN1);M(yN1,xN2),M(yN2,xN2)]; %Matrix mit den Werten der Stuetzpunkte
p = [(xN2-xN(i,j))/(xN2-xN1);(xN(i,j)-xN1)/(xN2-xN1)]; %Vektor mit Gewichtungsfunktionen fuer x-Richt.
o = [(yN2-yN(i,j))/(yN2-yN1),(yN(i,j)-yN1)/(yN2-yN1)]; %Vektor mit Gewichtungsfunktionen fuer y-Richt.
H1(round(yO(j)),round(xO(i))) = o*(Q*p); %Bilineare Interpolation in Matrixschreibweise
end end end end end toc
xN = bsxfun(@plus,x'*co,y*si)+origin(1);
yN = bsxfun(@minus,y*co,x'*si)+origin(2);
tic
H = zeros(240,320);
r_xO = round(xO);
r_yO = round(yO);
f_xN = floor(xN);
c_xN = ceil(xN);
f_yN = floor(yN);
c_yN = ceil(yN);
for i=1:lx
for j=1:ly
ifround(yO(j)) < sizeY && round(xO(i)) < sizeX && round(yO(j)) > 1 && round(xO(i)) > 1 if yN(i,j) < sizeY && xN(i,j) < sizeX && yN(i,j) > 1 && xN(i,j) > 1
xN1 = f_xN(i,j); %x-Koord. des1. Stuetzpunktes
xN2 = c_xN(i,j); %x-Koord. des2. Stuetzpunktes
yN1 = f_yN(i,j); %y-Koord. des1. Stuetzpunktes
yN2 = c_yN(i,j); %y-Koord. des2. Stuetzpunktes if xN1 == xN2 || yN1 == yN2
H(round(yO(j)),round(xO(i))) = M(yN1,xN1);
else
a1 = xN2-xN(i,j);
a2 = xN(i,j)-xN1;
b1 = yN2-yN(i,j);
b2 = yN(i,j)-yN1;
H(r_yO(j),r_xO(i)) = (b1*(M(yN1,xN1)*a1+M(yN2,xN1)*a2)+b2*(M(yN1,xN2)*a1+M(yN2,xN2)*a2))/((xN2-xN1)*(yN2-yN1));
end end end end end toc
xN = bsxfun(@plus,x'*co,y*si)+origin(1);
yN = bsxfun(@minus,y*co,x'*si)+origin(2);
tic
H = zeros(240,320);
f_xN = floor(xN);
c_xN = ceil(xN);
f_yN = floor(yN);
c_yN = ceil(yN);
r_xO = round(xO);
r_yO = round(yO);
match_O = r_yO < sizeY & r_yO > 1;
for i = 1:lx
if r_xO(i) < sizeX && r_xO(i) > 1
match = and(match_O,yN(i,:) < sizeY & xN(i,:) < sizeX & yN(i,:) > 1 & xN(i,:) > 1);
for j = 1:ly
if match(j)
xN1 = f_xN(i,j); %x-Koord. des1. Stuetzpunktes
xN2 = c_xN(i,j); %x-Koord. des2. Stuetzpunktes
yN1 = f_yN(i,j); %y-Koord. des1. Stuetzpunktes
yN2 = c_yN(i,j); %y-Koord. des2. Stuetzpunktes if xN1 == xN2 || yN1 == yN2
H(round(yO(j)),round(xO(i))) = M(yN1,xN1);
else
a1 = xN2-xN(i,j);
a2 = xN(i,j)-xN1;
b1 = yN2-yN(i,j);
b2 = yN(i,j)-yN1;
H(r_yO(j),r_xO(i)) = (b1*(M(yN1,xN1)*a1+M(yN2,xN1)*a2)+b2*(M(yN1,xN2)*a1+M(yN2,xN2)*a2))/((xN2-xN1)*(yN2-yN1));
end end end end end toc
Die erste Version benötigt bei mir ca. 2.68 s, die zweite ist 25 mal schneller und die dritte ist sogar 29 mal schneller als die erste Version (mit der oben im Code erstellten Matrix getestet).
Interessant ist, dass hier De-Vektorisieren der richtige Weg ist. Kannst du vielleicht ein paar Richtlinien nennen, wann Vektorisieren einer Schleife vorzuziehen ist, bzw. wann Vektorisieren einen zu hohen Datenverkehr oder Speicherkapazität erfordert?
Ich gebe mal hier den Hinweis, dass Code, der länger als eine Bildschirmseite ist, dringend (!) überarbeitet werden sollte, denn das kann niemand mehr nach 2 Jahren lesen bzw. verstehen. Vor allem können Kommentare helfen und ich rate ganz dringend dazu, dass man das einmal inhaltlich strukturiert
Während ich ansonsten jeden Sinnabschnitt mit einem Kommentar versehe und mit Leerzeilen abtrenne, find ich hier im Forum oft einen möglichst mageren Code sinnvoll. Ansonsten muss ich so viel scrollen.
Deshalb hatte ich mir auch die freiheit genommen, Kommentare aus Thomas' Orginalcode zu entfernen, um das eigentliche Problem möglichst kompakt behandeln zu können.
@Jan: Dein Code ist überschaubar, ich habe mich auf den Code im vorhergehenden Post von Thomas bezogen. Der Code von Thomas in der Forendarstellung ist mehr als eine Bildschirmseite bei Full HD Auflösung und da fehlen noch Kommentare usw
Es sind bei mir 2 B
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS