diff --git a/test_parallel_performance/PIVlab_commandline_parallel_speed.m b/test_parallel_performance/PIVlab_commandline_parallel_speed.m
new file mode 100644
index 0000000..198d101
--- /dev/null
+++ b/test_parallel_performance/PIVlab_commandline_parallel_speed.m
@@ -0,0 +1,231 @@
+function [calculation_times]=PIVlab_commandline_parallel_speed(directory,cpu_core_number)
+
+nr_of_cores =cpu_core_number; % integer, 1 means single core, greater than 1 means parallel
+if nr_of_cores > 1
+	try
+		local_cluster=parcluster('local'); % single node
+		corenum =  local_cluster.NumWorkers ; % fix : get the number of cores available
+	catch
+		warning('on');
+		warning('parallel local cluster can not be created, assigning number of cores to 1');
+		nr_of_cores = 1;
+	end
+end
+
+if nr_of_cores > 1
+	if pivparpool('size')~=nr_of_cores
+		pivparpool('close');
+		pivparpool('open',nr_of_cores);
+	end
+end
+
+%% Create list of images inside user specified directory
+
+% default directory: PIVlab/Examples
+addpath(fileparts(mfilename('fullpath')));
+%disp(fileparts(mfilename('fullpath')))
+
+suffix='*.jpg'; %*.bmp or *.tif or *.jpg or *.tiff or *.jpeg
+%disp(['Looking for ' suffix ' files in the selected directory.']);
+direc = dir([directory,filesep,suffix]); filenames={};
+[filenames{1:length(direc),1}] = deal(direc.name);
+filenames = sortrows(filenames); %sort all image files
+amount = length(filenames);
+
+%% Standard PIV Settings
+s = cell(15,2); % To make it more readable, let's create a "settings table"
+%Parameter                          %Setting           %Options
+s{1,1}= 'Int. area 1';              s{1,2}=128;         % window size of first pass
+s{2,1}= 'Step size 1';              s{2,2}=64;         % step of first pass
+s{3,1}= 'Subpix. finder';           s{3,2}=1;          % 1 = 3point Gauss, 2 = 2D Gauss
+s{4,1}= 'Mask';                     s{4,2}=[];         % If needed, generate via: imagesc(image); [temp,Mask{1,1},Mask{1,2}]=roipoly;
+s{5,1}= 'ROI';                      s{5,2}=[];         % Region of interest: [x,y,width,height] in pixels, may be left empty
+s{6,1}= 'Nr. of passes';            s{6,2}=2;          % 1-4 nr. of passes
+s{7,1}= 'Int. area 2';              s{7,2}=64;         % second pass window size
+s{8,1}= 'Int. area 3';              s{8,2}=16;         % third pass window size
+s{9,1}= 'Int. area 4';              s{9,2}=16;         % fourth pass window size
+s{10,1}='Window deformation';       s{10,2}='*linear'; % '*spline' is more accurate, but slower
+s{11,1}='Repeated Correlation';     s{11,2}=0;         % 0 or 1 : Repeat the correlation four times and multiply the correlation matrices.
+s{12,1}='Disable Autocorrelation';  s{12,2}=0;         % 0 or 1 : Disable Autocorrelation in the first pass.
+s{13,1}='Correlation style';        s{13,2}=0;         % 0 or 1 : Use circular correlation (0) or linear correlation (1).
+s{14,1}='Repeat last pass';   s{14,2}=0; % 0 or 1 : Repeat the last pass of a multipass analyis
+s{15,1}='Last pass quality slope';   s{15,2}=0.025; % Repetitions of last pass will stop when the average difference to the previous pass is less than this number.
+
+
+%% Standard image preprocessing settings
+p = cell(10,1);
+%Parameter                       %Setting           %Options
+p{1,1}= 'ROI';                   p{1,2}=s{5,2};     % same as in PIV settings
+p{2,1}= 'CLAHE';                 p{2,2}=1;          % 1 = enable CLAHE (contrast enhancement), 0 = disable
+p{3,1}= 'CLAHE size';            p{3,2}=50;         % CLAHE window size
+p{4,1}= 'Highpass';              p{4,2}=0;          % 1 = enable highpass, 0 = disable
+p{5,1}= 'Highpass size';         p{5,2}=15;         % highpass size
+p{6,1}= 'Clipping';              p{6,2}=0;          % 1 = enable clipping, 0 = disable
+p{7,1}= 'Wiener';                p{7,2}=0;          % 1 = enable Wiener2 adaptive denoise filter, 0 = disable
+p{8,1}= 'Wiener size';           p{8,2}=3;          % Wiener2 window size
+p{9,1}= 'Minimum intensity';     p{9,2}=0.0;        % Minimum intensity of input image (0 = no change)
+p{10,1}='Maximum intensity';     p{10,2}=1.0;       % Maximum intensity on input image (1 = no change)
+
+%% other settings
+pairwise = 1; % 0 for [A+B], [B+C], [C+D]... sequencing style, and 1 for [A+B], [C+D], [E+F]... sequencing style
+
+%% PIV analysis loop
+if pairwise == 1
+	if mod(amount,2) == 1 %Uneven number of images?
+		disp('Image folder should contain an even number of images.')
+		%remove last image from list
+		amount=amount-1;
+		filenames(size(filenames,1))=[];
+	end
+	
+	%disp(['Found ' num2str(amount) ' images (' num2str(amount/2) ' image pairs).'])
+	x=cell(amount/2,1);
+	y=x;
+	u=x;
+	v=x;
+else
+	%disp(['Found ' num2str(amount) ' images'])
+	x=cell(amount-1,1);
+	y=x;
+	u=x;
+	v=x;
+end
+
+typevector=x; %typevector will be 1 for regular vectors, 0 for masked areas
+correlation_map=x; % correlation coefficient
+
+%% Pre-load the image names out side of the parallelizable loop
+slicedfilename1=cell(0);
+slicedfilename2=cell(0);
+j = 1;
+for i=1:1+pairwise:amount-1
+	slicedfilename1{j}=filenames{i}; % begin
+	slicedfilename2{j}=filenames{i+1}; % end
+	j = j+1;
+end
+
+dummy=imread(fullfile(directory,filenames{1}));
+disp(['    Current image size:    ' num2str(size(dummy,1)) ' px'])
+
+
+tstart=tic;
+%% Main PIV analysis loop:
+% parallel
+disp('Please wait....')
+if nr_of_cores > 1
+	
+	
+	
+	parfor i=1:size(slicedfilename1,2)  % index must increment by 1
+		
+		[x{i}, y{i}, u{i}, v{i}, typevector{i},correlation_map{i}] = ...
+			piv_analysis(directory, slicedfilename1{i}, slicedfilename2{i},p,s,nr_of_cores,false);
+	end
+else % sequential loop
+	
+	for i=1:size(slicedfilename1,2)  % index must increment by 1
+		
+		[x{i}, y{i}, u{i}, v{i}, typevector{i},correlation_map{i}] = ...
+			piv_analysis(directory, slicedfilename1{i}, slicedfilename2{i},p,s,nr_of_cores,false);
+		
+		%disp([int2str((i+1)/amount*100) ' %']);
+		
+	end
+end
+
+
+%% PIV postprocessing loop
+% Standard image post processing settings
+
+r = cell(6,1);
+%Parameter     %Setting                                     %Options
+r{1,1}= 'Calibration factor, 1 for uncalibrated data';      r{1,2}=1;                   % Calibration factor for u
+r{2,1}= 'Calibration factor, 1 for uncalibrated data';      r{2,2}=1;                   % Calibration factor for v
+r{3,1}= 'Valid velocities [u_min; u_max; v_min; v_max]';    r{3,2}=[-50; 50; -50; 50];  % Maximum allowed velocities, for uncalibrated data: maximum displacement in pixels
+r{4,1}= 'Stdev check?';                                     r{4,2}=1;                   % 1 = enable global standard deviation test
+r{5,1}= 'Stdev threshold';                                  r{5,2}=7;                   % Threshold for the stdev test
+r{6,1}= 'Local median check?';                              r{6,2}=1;                   % 1 = enable local median test
+r{7,1}= 'Local median threshold';                           r{7,2}=3;                   % Threshold for the local median test
+
+u_filt=cell(size(u));
+v_filt=cell(size(v));
+typevector_filt=typevector;
+
+if nr_of_cores >1 % parallel
+	
+	
+	if pivparpool('size')~=nr_of_cores
+		pivparpool('open',nr_of_cores);
+	end
+	
+	
+	parfor PIVresult=1:size(x,1)
+		
+		[u_filt{PIVresult,1}, v_filt{PIVresult,1},typevector_filt{PIVresult,1}]= ...
+			post_proc_wrapper(u{PIVresult,1},v{PIVresult,1},typevector{PIVresult,1},r,true);
+		
+	end
+	
+else % sequential loop
+	
+	for PIVresult=1:size(x,1)
+		
+		[u_filt{PIVresult,1}, v_filt{PIVresult,1},typevector_filt{PIVresult,1}]= ...
+			post_proc_wrapper(u{PIVresult,1},v{PIVresult,1},typevector{PIVresult,1},r,true);
+		
+	end
+	
+end
+calculation_times=toc(tstart);
+%% clean up parallel pool, and cluster
+%{
+if nr_of_cores >1 % parallel
+	poolobj = gcp('nocreate'); % GET the current parallel pool
+	if ~isempty(poolobj ); delete(poolobj );end
+	clear local_cluster;
+end
+%}
+%%
+%save(fullfile(directory, [filenames{1} '_' filenames{end} '_' num2str(amount) '_frames_result_.mat']));
+
+%%
+%clearvars -except p s r x y u v typevector directory filenames u_filt v_filt typevector_filt correlation_map
+disp('DONE.')
+
+end
+function [u_filt, v_filt,typevector_filt] = post_proc_wrapper(u,v,typevector,post_proc_setting,paint_nan)
+% wrapper function for PIVlab_postproc
+
+% INPUT
+% u, v: u and v components of vector fields
+% typevector: type vector
+% post_proc_setting: post processing setting
+% paint_nan: bool, whether to interpolate missing data
+
+% OUTPUT
+% u_filt, v_filt: post-processed u and v components of vector fields
+% typevector_filt: post-processed type vector
+
+
+[u_filt,v_filt] = PIVlab_postproc(u,v, ...
+	post_proc_setting{1,2},...
+	post_proc_setting{2,2},...
+	post_proc_setting{3,2},...
+	post_proc_setting{4,2},...
+	post_proc_setting{5,2},...
+	post_proc_setting{6,2},...
+	post_proc_setting{7,2});
+
+typevector_filt = typevector; % initiate
+typevector_filt(isnan(u_filt))=2;
+typevector_filt(isnan(v_filt))=2;
+typevector_filt(typevector==0)=0; %restores typevector for mask
+
+if paint_nan
+	u_filt=inpaint_nans(u_filt,4);
+	v_filt=inpaint_nans(v_filt,4);
+end
+
+
+end
+
diff --git a/test_parallel_performance/Williams_results.jpg b/test_parallel_performance/Williams_results.jpg
new file mode 100644
index 0000000..15caf99
Binary files /dev/null and b/test_parallel_performance/Williams_results.jpg differ
diff --git a/test_parallel_performance/do_parallel_performance_analysis.m b/test_parallel_performance/do_parallel_performance_analysis.m
new file mode 100644
index 0000000..c96d7a5
--- /dev/null
+++ b/test_parallel_performance/do_parallel_performance_analysis.m
@@ -0,0 +1,40 @@
+%Parallel computing performance analysis including post processing.
+
+%% Add PIVlab path so we can access functions
+addpath('C:\Users\trash\Documents\MATLAB\PIVlab')
+
+%% Settings
+image_sizes=[600 1200 2500];
+cores_to_use=[1 3  5 7 9];
+nr_of_image_pairs=50;
+
+%% 1: generate piv images with different sizes
+for image_size=image_sizes
+	disp(['         Image size:         ' num2str(image_size) ' px'])
+	mkdir (num2str(image_size));
+	generate_images(fullfile(pwd,num2str(image_size)),nr_of_image_pairs,image_size)
+end
+
+%% 2: analyze these images with different amount of cores
+idx2=1;
+calculation_times=zeros(numel(image_sizes),numel(cores_to_use));
+for cores=cores_to_use
+	idx=1;
+	for image_size=image_sizes
+		calculation_times(idx,idx2)=PIVlab_commandline_parallel_speed(fullfile(pwd,num2str(image_size)),cores);
+		idx=idx+1;
+	end
+	idx2=idx2+1;
+end
+
+%% 3: Plot the results
+%rows=image size
+%cols=num of cores
+figure;
+semilogy(calculation_times','linewidth',2)
+legend(compose('%d px',image_sizes))
+grid on
+xlabel('nr of cores')
+set(gca,'xticklabel',compose('%d cores',cores_to_use))
+ylabel('processing time')
+
diff --git a/test_parallel_performance/generate_images.m b/test_parallel_performance/generate_images.m
new file mode 100644
index 0000000..97cf018
--- /dev/null
+++ b/test_parallel_performance/generate_images.m
@@ -0,0 +1,136 @@
+function generate_images (outputpath,amount,size)
+selpath = outputpath;
+noise = 0.0005; %noise 0 bis 0.05
+z_move=10; %z_move
+partAm=size^2*0.01;
+Z=0.333; %sheet thickness
+dt=3; %particle diameter
+ddt=1; %particle diameter variation
+
+for iiii = 0:amount-1
+disp(['Progress: ' int2str(iiii/amount*100) ' %']);
+	%% Generate random artificial particle images
+	u = rand(1)*ones(size,size);
+	v = rand(1)*ones(size,size);
+	offsety=v; %zero displacement in y direction
+	offsetx=u; %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=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
+