Shrediquette / PIVlab

Particle Image Velocimetry for Matlab, official repository
https://shrediquette.github.io/PIVlab/
MIT License
126 stars 30 forks source link

how to generate 10000 synthetic pics? #44

Closed WangJieHHU closed 3 years ago

WangJieHHU commented 3 years ago

If i want to generate 10000 pairs of synthetic pictures,which func should i use?

Shrediquette commented 3 years ago

This will wok:


clear all; clc;drawnow selpath = uigetdir; noise = 0; %noise 0 bis 0.05 z_move=0; %z_move size=600; %image size x and y partAm=20000; Z=0.333; %sheet thickness dt=4; %particle diameter ddt=2; %particle diameter variation disp(['Generating random artificial PIV images with ' num2str(partAm) ' particles...']) for iiii = 0:99 disp(['Progress: ' int2str(iiii/99100) ' %']); %% Generate random artificial particle images u = zeros(size,size); v=zeros(size,size); offsety=v+0; %zero displacement in y direction offsetx=u+iiii/10; %uniform x displacement increases from 0 to 99/10 pixels [x,y]=meshgrid(1:1:size); i=[]; j=[]; sizey=size; sizex=size; A=zeros(sizey,sizex); B=A; z0_pre=randn(partAm,1); %normal distributed sheet intensity randn('state', sum(100clock)); %#ok<RAND> z1_pre=randn(partAm,1); %normal distributed sheet intensity z0=z0_pre(z_move/200+0.5)+z1_pre(1-((z_move/200+0.5))); z1=z1_pre(z_move/200+0.5)+z0_pre*(1-((z_move/200+0.5)));

I0=255*exp(-(Z^2./(0.125*z0.^2))); %particle intensity
I0(I0>255)=255;
I0(I0<0)=0;

I1=255*exp(-(Z^2./(0.125*z1.^2))); %particle intensity
I1(I1>255)=255;
I1(I1<0)=0;

randn('state', sum(100*clock));
d=randn(partAm,1)/2; %particle diameter distribution
d=dt+d*ddt;
d(d<0)=0;
rand('state', sum(100*clock));
x0=rand(partAm,1)*sizex;
y0=rand(partAm,1)*sizey;
rd = -8.0 ./ d.^2;

xlimit1=floor(x0-d/2); %x min particle extent image1
xlimit2=ceil(x0+d/2); %x max particle extent image1
ylimit1=floor(y0-d/2); %y min particle extent image1
ylimit2=ceil(y0+d/2); %y max particle extent image1
xlimit2(xlimit2>sizex)=sizex;
xlimit1(xlimit1<1)=1;
ylimit2(ylimit2>sizey)=sizey;
ylimit1(ylimit1<1)=1;

%calculate particle extents for image2 (shifted image)
x0integer=round(x0);
x0integer(x0integer>sizex)=sizex;
x0integer(x0integer<1)=1;
y0integer=round(y0);
y0integer(y0integer>sizey)=sizey;
y0integer(y0integer<1)=1;

xlimit3=zeros(partAm,1);
xlimit4=xlimit3;
ylimit3=xlimit3;
ylimit4=xlimit3;
for n=1:partAm
    xlimit3(n,1)=floor(x0(n)-d(n)/2-offsetx((y0integer(n)),(x0integer(n)))); %x min particle extent image2
    xlimit4(n,1)=ceil(x0(n)+d(n)/2-offsetx((y0integer(n)),(x0integer(n)))); %x max particle extent image2
    ylimit3(n,1)=floor(y0(n)-d(n)/2-offsety((y0integer(n)),(x0integer(n)))); %y min particle extent image2
    ylimit4(n,1)=ceil(y0(n)+d(n)/2-offsety((y0integer(n)),(x0integer(n)))); %y max particle extent image2
end
xlimit3(xlimit3<1)=1;
xlimit4(xlimit4>sizex)=sizex;
ylimit3(ylimit3<1)=1;
ylimit4(ylimit4>sizey)=sizey;

ctr=0;
for n=1:partAm
    ctr=ctr+1;
    if ctr==10000
        ctr=0;
        %fprintf('.')
    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)+I1(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

%% Create random Background "glow"
bg=im2bw(rand(size,size),0.99999);
SE = strel('line',round(rand*90)+90,round(rand*180));
bg2=imdilate(bg,SE);
bg3 = imgaussfilt(double(bg2),40);
bg3=bg3/max(max(bg3));
bg3=bg3*0.1;
bg=im2bw(rand(size,size),0.999997);
SE = strel('disk',round(rand*100)+200,4);
bg2=imdilate(bg,SE);
bg4 = 0.12*imgaussfilt(double(bg2),100);
bg3(isnan(bg3))=0;
bg4(isnan(bg4))=0;

A=A+bg3*rand*255+bg4*rand*255;
B=B+bg3*rand*255+bg4*rand*255;

A(A>255)=255;
B(B>255)=255;
A=imnoise(uint8(A),'gaussian',0,noise);
B=imnoise(uint8(B),'gaussian',0,noise);

A=uint8(A);
B=uint8(B);

%True velocities, save if desired
x_real=x;
y_real=y;
u_real=u;
v_real=v;

imwrite(A,[selpath '\synth_image_' sprintf('%5.5d',iiii) '_A.jpg'],'Quality',97);
imwrite(B,[selpath '\synth_image_' sprintf('%5.5d',iiii) '_B.jpg'],'Quality',97);

end

Shrediquette commented 3 years ago

... but code formatting is broken in this github website

WangJieHHU commented 3 years ago

This will work: `clear all; clc;drawnow selpath = uigetdir; noise = 0; %noise 0 bis 0.05 z_move=0; %z_move size=600; %image size x and y partAm=20000; Z=0.333; %sheet thickness dt=4; %particle diameter ddt=2; %particle diameter variation disp(['Generating random artificial PIV images with ' num2str(partAm) ' particles...']) for iiii = 0:99 disp(['Progress: ' int2str(iiii/99_100) ' %']); %% Generate random artificial particle images u = zeros(size,size); v=zeros(size,size); offsety=v+0; %zero displacement in y direction offsetx=u+iiii/10; %uniform x displacement increases from 0 to 99/10 pixels [x,y]=meshgrid(1:1:size); i=[]; j=[]; sizey=size; sizex=size; A=zeros(sizey,sizex); B=A; z0_pre=randn(partAm,1); %normal distributed sheet intensity randn('state', sum(100_clock)); %#ok<_RAND> z1_pre=randn(partAm,1); %normal distributed sheet intensity z0=z0pre(z_move/200+0.5)+z1_pre(1-((z_move/200+0.5))); z1=z1_pre(z_move/200+0.5)+z0_pre*(1-((z_move/200+0.5)));

I0=255*exp(-(Z^2./(0.125*z0.^2))); %particle intensity
I0(I0>255)=255;
I0(I0<0)=0;

I1=255*exp(-(Z^2./(0.125*z1.^2))); %particle intensity
I1(I1>255)=255;
I1(I1<0)=0;

randn('state', sum(100*clock));
d=randn(partAm,1)/2; %particle diameter distribution
d=dt+d*ddt;
d(d<0)=0;
rand('state', sum(100*clock));
x0=rand(partAm,1)*sizex;
y0=rand(partAm,1)*sizey;
rd = -8.0 ./ d.^2;

xlimit1=floor(x0-d/2); %x min particle extent image1
xlimit2=ceil(x0+d/2); %x max particle extent image1
ylimit1=floor(y0-d/2); %y min particle extent image1
ylimit2=ceil(y0+d/2); %y max particle extent image1
xlimit2(xlimit2>sizex)=sizex;
xlimit1(xlimit1<1)=1;
ylimit2(ylimit2>sizey)=sizey;
ylimit1(ylimit1<1)=1;

%calculate particle extents for image2 (shifted image)
x0integer=round(x0);
x0integer(x0integer>sizex)=sizex;
x0integer(x0integer<1)=1;
y0integer=round(y0);
y0integer(y0integer>sizey)=sizey;
y0integer(y0integer<1)=1;

xlimit3=zeros(partAm,1);
xlimit4=xlimit3;
ylimit3=xlimit3;
ylimit4=xlimit3;
for n=1:partAm
  xlimit3(n,1)=floor(x0(n)-d(n)/2-offsetx((y0integer(n)),(x0integer(n)))); %x min particle extent image2
  xlimit4(n,1)=ceil(x0(n)+d(n)/2-offsetx((y0integer(n)),(x0integer(n)))); %x max particle extent image2
  ylimit3(n,1)=floor(y0(n)-d(n)/2-offsety((y0integer(n)),(x0integer(n)))); %y min particle extent image2
  ylimit4(n,1)=ceil(y0(n)+d(n)/2-offsety((y0integer(n)),(x0integer(n)))); %y max particle extent image2
end
xlimit3(xlimit3<1)=1;
xlimit4(xlimit4>sizex)=sizex;
ylimit3(ylimit3<1)=1;
ylimit4(ylimit4>sizey)=sizey;

ctr=0;
for n=1:partAm
  ctr=ctr+1;
  if ctr==10000
      ctr=0;
      %fprintf('.')
  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)+I1(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

%% Create random Background "glow"
bg=im2bw(rand(size,size),0.99999);
SE = strel('line',round(rand*90)+90,round(rand*180));
bg2=imdilate(bg,SE);
bg3 = imgaussfilt(double(bg2),40);
bg3=bg3/max(max(bg3));
bg3=bg3*0.1;
bg=im2bw(rand(size,size),0.999997);
SE = strel('disk',round(rand*100)+200,4);
bg2=imdilate(bg,SE);
bg4 = 0.12*imgaussfilt(double(bg2),100);
bg3(isnan(bg3))=0;
bg4(isnan(bg4))=0;

A=A+bg3*rand*255+bg4*rand*255;
B=B+bg3*rand*255+bg4*rand*255;

A(A>255)=255;
B(B>255)=255;
A=imnoise(uint8(A),'gaussian',0,noise);
B=imnoise(uint8(B),'gaussian',0,noise);

A=uint8(A);
B=uint8(B);

%True velocities, save if desired
x_real=x;
y_real=y;
u_real=u;
v_real=v;

imwrite(A,[selpath '\synth_image_' sprintf('%5.5d',iiii) '_A.jpg'],'Quality',97);
imwrite(B,[selpath '\synth_image_' sprintf('%5.5d',iiii) '_B.jpg'],'Quality',97);

end `

Thank u!I‘m going to try it!