Verfasst am: 11.11.2010, 12:13
Titel: For loop vektorisieren...
Hi!
Ich möchte zwei synthetische Partikelbilder (A & B) mit Matlab erstellen. Nur leider dauert das erstellen bei ganz kleinen Testbildern (200x200 px) mit nur 1000 Partikeln bereits 4 Sekunden. Am Ende möchte ich gerne 1000 x 1000 px große Bilder mit 15000 Partikeln erstellen. Das dürfte ja mehrere Tage dauern...
Folgenden Code habe ich geschrieben:
Code:
A=zeros(200,200);
B=A;
partAm=1000; % Amount of Particles
Z=0.2; % sheet thickness
dt=20; % particle diameter
ddt=10; % particle diameter variation
z0=abs(randn(partAm,1)); % normal distributed sheet intensity
I0=255*exp(-(Z^2./(0.125*z0))); % particle intensity
I0(I0>255)=255;
I0(I0<0)=0;
d=randn(partAm,1)/2; % particle diameter distribution
d=dt+d*ddt;
d(d<0)=0;
x0=rand(partAm,1)*size(A,2); % x position of the particles in A
y0=rand(partAm,1)*size(A,1); % y position of the particles in A
offsetx=0.2; % x offset for image B
offsety=2; % y offset for image B tic for n=1:partAm % for each particle for i=1:size(A,1) for j=1:size(A,2)
A(i,j)=A(i,j)+I0(n)*exp((-(i-x0(n))^2-(j-y0(n))^2)/((1/8*d(n)))); % place particle with gaussian intensity profile
B(i,j)=B(i,j)+I0(n)*exp((-(i-x0(n)+offsetx)^2-(j-y0(n)+offsety)^2)/((1/8*d(n))));
if A(i,j)>255
A(i,j)=255;
end if B(i,j)>255
B(i,j)=255;
end end end end toc imshow(uint8(A)) figure;imshow(uint8(B))
Probiert ihn doch einfach mal aus um zu sehen was ich mit "Partikelbilder" meine.
Wichtig ist hier nur der for loop (besonders die inneren beiden Schleifen). Für jeden der kleinen Partikel wird jedesmal das ganze Bild durchlaufen. Das dauert natürlich ewig. Fällt jemandem von euch ein wie man das schneller machen könnte? Kann das nicht irgendwie vektorisiert werden...?
Danke schonmal für eure Tipps & viele Grüße,
William
_________________
Verfasst am: 11.11.2010, 14:34
Titel: Re: For loop vektorisieren...
Hallo Willa,
Ich empfehle vor der Vektorisierung erstmal die Standard-Methoden der effizienten Programmierung anzuwenden:
1. So viele Berechnungen aus den Schleifen auslagern wie möglich.
In Deinem Beipiel wird "1/8*d(n)" in jeder Iteration neu berechnet, dabei ändert sich der Wert ja nur in der äußersten Schleife. Also definiere eine temporäre Variable dort:
2. Die Multiplikation ist schneller als die Division. Also definiere gleich den reziproken Wert der temporären Variable:
Code:
... for n=1:partAm
r = 8.0 / d(n); % Oder diese Rechnung auch gleich nach aussen! ...
A(i,j)=A(i,j)+I0(n)*exp((-(i-x0(n))^2-(j-y0(n))^2) * r);
3. Wenn man ein Array elementweise bearbeitet, kann auf die Elemente schneller zugegriffen werden, die im Speicher benachbart sind, also die erste Dimension. Deshalb wird die Schleife über die erste Dimension nach innen gelegt, und die über die zweite Dimension nach aussen. Also "for i" und "for j" vertauschen.
4. Der Ausdruck "(j-y0(n))^2)" ändert sich danach auch nicht mehr in der inneren Schleife, also wird das nur noch in der "for j"-Schleife berechnet und wieder temporär gespeichert.
5. SIZE ist ein Befehl, der Zeit benötigt. Es ist daher schneller die Dimension der Matrix einmal abzufragen und temporär zu speichern.
6. Die Begrenzung der maximalen Werte auf 255 kann auch einmalig zum Schluß durchgeführt werden.
s1 = 200;
s2 = 200;
A = zeros(s1,s2);
B = zeros(s1,s2);
partAm=1000; % Amount of Particles
Z=0.2; % sheet thickness
dt=20; % particle diameter
ddt=10; % particle diameter variation
z0=abs(randn(partAm,1)); % normal distributed sheet intensity
I0=255*exp(-(Z^2./(0.125*z0))); % particle intensity
I0(I0>255)=255;
I0(I0<0)=0;
d=randn(partAm,1)/2; % particle diameter distribution
d=dt+d*ddt;
d(d<0)=0;
x0=rand(partAm,1)*s2; % x position of the particles in A
y0=rand(partAm,1)*s1; % y position of the particles in A
offsetx=0.2; % x offset for image B
offsety=2; % y offset for image B
rd = -8.0 ./ d;
tic for n=1:partAm % for each particle
r = rd(n);
for j = 1:s2
rj = (j-y0(n))^2;
rjo = (j-y0(n)+offsety)^2;
for i = 1:s1
ri = -(i-x0(n))^2; % Nur zur Übersichtlichkeit temporär
rio = -(i-x0(n)+offsetx)^2;
A(i,j) = A(i,j) + I0(n) * exp((ri - rj) * r);
B(i,j) = B(i,j) + I0(n) * exp((rio -rjo) * r);
end end end
A(A > 255) = 255;
B(B > 255) = 255;
toc imshow(uint8(A)) figure;imshow(uint8(B))
Damit kommt ich auf meinem betagten 1.5GHz PentiumM schon mal von 9.8sec auf 5.7sec. Das ist zwar nicht umwerfend viel schneller, aber für das einfache Umsortieren von Befehlen schon ganz nett. Zudem ist der Code (meiner Meinung nach) übersichtlicher geworden, was das Vektorisieren erleichtert.
Die Multiplikation mit "I0(n)" nach den Schleifen auszuführen bringt übrigens nichts. Matlab's JIT-Compiler ist da bereits sehr effizient.
Eine Vektorisierung schicke ich im nächsten Posting.
Hi Jan,
vielen Dank für deine ausführliche Antwort! Ich werde deine Vorschläge schnellstmöglich umsetzen. Mittlerweile habe ich einen Weg gefunden, der meinen Code 173 mal schneller laufen lässt. Ich gehe nicht das ganze Bild durch für jeden Partikel, sondern immer nur einen kleinen Bereich (setzt sich aus Partikelgröße und Versatz zusammen) um den zu zeichnenden Partikel. mal sehen wie schnell es wird wenn ich deine Vorschläge umgesetzt habe....
_________________
Dies benötigt nun 4.8sec, statt 5.7 für die Version mit der doppelten Schleife. Offenbar kann Matlab die Schleifen schon ziemlich gut abarbeiten! Wo ist jetzt der Flaschenhals? In EXP!
Die EXP-Funktion wird für eine Matrix ausgerechnet, aber man kann sie separieren und auf die einzelnen Vektoren anwenden. Dann muss man nur noch 200+200 EXP-Werte berechnen, und nicht 200*200!
2.1 sec! (9.8 in der Original-Version, 5.7 mit effizienten Schleifen, 4.8 mit BSXFUN statt der Schleifen).
Faktor 4.5 ist zwar ganz nett, aber der Übergang von 200x200 auf 1000x1000 Matrizen wird dadurch nicht gerade lustig. Wie könnte man das Programm dann wirklich effizient machen?
Die Summen A+I0(n)*X addieren hauptsächlich Nullen zu Nullen! Dann wäre es viel besser, man würde die maximale Ausdehnung der "Punkte" zunerst ausrechnen. Wenn sie maximal 100 Pixel Ausdehnung haben, muss man die neuen Werte ja auch nur zu den entsprechenden Sub-Matrix addieren, und auch weit weniger EXP-Werte berechnen. Bei den 200x200 Matrizen bringt das wohl einen Faktor 4, bei 1000X1000 aber schon ein Faktor 100 (5*5*4)!
Jetzt läuft der Code 388 mal schneller als zu Beginn, ich habe versucht die Kommentare von Jan zu berücksichtigen. In meiner (bzw. "der Hochschule ihm seiner") Version (7.1 R14) gibts leider noch kein bsxfun. Den größten Geschwindigkeitszuwachs bringt wie erwartet die Begrenzung bei der Errechnung der Partikelintensitäten.
Für ein Bild der Größe 1000 * 1000 px bei 100000 Partikeln braucht mein Computer nur noch 1.21 Sekunden
Hier der Code:
for n=1:partAm
xlimit3(n,1)=floor(x0(n)-d(n)/2-offsetx(round(y0integer(n)),round(x0integer(n)))); %x min particle extent image2
xlimit4(n,1)=ceil(x0(n)+d(n)/2-offsetx(round(y0integer(n)),round(x0integer(n)))); %x max particle extent image2
ylimit3(n,1)=floor(y0(n)-d(n)/2-offsety(round(y0integer(n)),round(x0integer(n)))); %y min particle extent image2
ylimit4(n,1)=ceil(y0(n)+d(n)/2-offsety(round(y0integer(n)),round(x0integer(n)))); %y max particle extent image2
if xlimit3(n)<1
xlimit3(n)=1;
end if xlimit4(n)>sizex
xlimit4(n)=sizex;
end if ylimit3(n)<1
ylimit3(n)=1;
end if ylimit4(n)>sizey
ylimit4(n)=sizey;
end
r = rd(n);
for j=xlimit1(n):xlimit2(n)
rj = (j-x0(n))^2;
for i=ylimit1(n):ylimit2(n)
A(i,j)=A(i,j)+I0(n)*exp((rj+(i-y0(n))^2)*r);
end end for j=xlimit3(n):xlimit4(n) for i=ylimit3(n):ylimit4(n)
B(i,j)=B(i,j)+I0(n)*exp((-(j-x0(n)+offsetx(i,j))^2-(i-y0(n)+offsety(i,j))^2)*-rd(n)); %place particle with gaussian intensity profile end end end toc
Sicherlich ist das noch deutlich weiter optimierbar (z.B. mit aktuelleren Matlab Versionen...), aber der Code läuft jetzt für meine Anwendung schnell genug.
Vielen dank an Jan & viele Grüße,
William
_________________
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