Verfasst am: 02.07.2008, 23:34
Titel: Künstliches Partikelbild für PIV erstellen
Hallo!
Ich habe kürzlich meine eigene Software für PIV (Particle Image Velocimetry) in Matlab programmiert. http://de.wikipedia.org/wiki/Particle_Image_Velocimetry Die funktioniert sehr schön, aber jetzt möchte ich mal überprüfen wie präzise das so ist.
Dafür, so dachte ich mir, erstelle ich mal 2 Bilder mit einem bekannten Partikelversatz. Da meine PIV-Software den Partikelversatz bis auf Subpixel-genauigkeit berechnet (per COG-weighted-to-gray Peakfinder), brauche ich 2 Bilder mit genau definiertem Subpixel-Versatz.
Ich habe folgende Funktion geschrieben:
Code:
function[A,B] = CreateParticleImage () % create a BIG image with BIG particles... (LargeImageSize) % Shift these particles an integer number of pixels % Scale down the Image (FinalImageSize) % ==> The particles will now be shifted a fraction of a pixel (Subpixel) % Do this whole thing for 10 "Layers" with different intensities % (This simulates the intensity profile of a laser sheet) clc;clear;
tic
LargeImageSize=2000;
FinalImageSize=400;
amountLayers=10; % 10 different intensities
A=zeros(LargeImageSize,LargeImageSize,amountLayers); % preallocate mem
for i = 1:amountLayers
rand('twister', sum(100*clock)); % shuffle the pot...
A(:,:,i)=rand(LargeImageSize,LargeImageSize); % random numbers Which=A(:,:,i)>0.9998; % Only pixels above this thresh will be kept
A(:,:,i)=Which.*A(:,:,i)*255;
H = fspecial('disk',7); % i know that this should be outside the loop...
A(:,:,i) = imfilter(A(:,:,i),H,'symmetric'); % make a "particle" from a single pixel
H = fspecial('average',7); % i know that this should be outside the loop...
A(:,:,i) = imfilter(A(:,:,i),H,'symmetric'); % simulate "scattering of the laser light"
A(:,:,i)=A(:,:,i)*i*255/27; % Layer1 -> low intensity, ... , Layer10 -> high intensity end for i =2:amountLayers
A(:,:,1)=A(:,:,1)+A(:,:,i); % combine the Layers end
A(:,:,2:amountLayers)=[]; % remove the rest of the layers
shiftx=1; % Shift the particles in X direction
shifty=13; % Shift the particles in Y direction
B = circshift(A, [shifty, shiftx]); % Do the displacement disp(['Shift in X: ' num2str(shiftx/(LargeImageSize/FinalImageSize)) ' px']); % how many subpixels was it shifted disp(['Shift in Y: ' num2str(shifty/(LargeImageSize/FinalImageSize)) ' px']); % how many subpixels was it shifted
A=imresize(A,[FinalImageSize FinalImageSize],'bicubic'); % scale down the image for subpixel shift
B=imresize(B,[FinalImageSize FinalImageSize],'bicubic'); % scale down the image for subpixel shift
A=A/(max(max(A)))*255; % normalize image
B=B/(max(max(B)))*255; % normalize image % imwrite(uint8(A),'artificial1.tif','tiff','Compression', 'none'); % imwrite(uint8(B),'artificial2.tif','tiff','Compression', 'none'); figure;imshow(uint8(A));
figure;imshow(uint8(B));
disp(['Wow, this took "only" ' num2str(toc) ' seconds....! :-(']) end
Ich habe versucht das ganze möglichst verständlich zu machen durch viele kommentare.
Das Problem an meinem Verfahren ist nun: Es dauert viiiieel zu lange (~36 sec. kann ich ja auch verstehen bei meiner Methode). Ich muss ganz viele Testbilder machen, da brauche ich ja Wochen.... Außerdem errechnet meine PIV Software aus den Bildern meist nicht den Wert den ich vorgebe. Ich vertraue aber meiner PIV Sofware mehr als meiner Function oben. Kann es sein dass ich den falschen Weg gewählt habe um "Subpixelverschobene" Partikelbilder zu generieren...?
Habt Ihr evtl. ein Paar Tipps was ich besser bzw. anders machen könnte? Wie ihr evtl. am Quelltext merkt bin ich noch nicht so der Profi... Aber ich versuche ständig mich zu verbessern.
Vielen Dank und viele Grüße,
William
Hi! Ich habs geschafft mein Problem selber zu lösen..... Statt quadratische einzelne Pixel zu erstellen und die danach zu blurren erstelle ich nun Normalverteile "Partikel" der Größe PartSize.
Es dauert nun nur noch 4 sekunden um 9 ebenen mit untersch. Intensität zu erstellen.
Code:
PartSize=7;
PartAm=500;
A=zeros(500,500);
i=1;
for j= 1:PartAm
row=ceil(rand*(size(A,1)-(PartSize-1)))+(PartSize-1)/2;
col=ceil(rand*(size(A,2)-(PartSize-1)))+(PartSize-1)/2;
for xdir=-(PartSize-1)/2:(PartSize-1)/2 for ydir =-(PartSize-1)/2:(PartSize-1)/2
x0=col;
y0=row;
x=col+xdir;
y=row+ydir;
A(y,x,i)=A(y,x,i)+255*exp((-(x-x0)^2-(y-y0)^2) / ((1/8)*PartSize^2)); % gaussian intensity particle profile end end end imshow(uint8(A))
Jetzt muss ich es nur noch schaffen in A alle punkte >255 auf 255 zu setzen...
Viele Grüße,
William
Einstellungen und Berechtigungen
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