From 6ea375a98ca5f02e6dd68a3c6749c93c8030f2fc Mon Sep 17 00:00:00 2001 From: William Thielicke <42935970+Shrediquette@users.noreply.github.com> Date: Thu, 19 Aug 2021 00:10:10 +0200 Subject: [PATCH] Add files via upload --- .../PIVlab_commandline_parallel_speed.m | 231 ++++++++++++++++++ .../Williams_results.jpg | Bin 0 -> 40541 bytes .../do_parallel_performance_analysis.m | 40 +++ test_parallel_performance/generate_images.m | 136 +++++++++++ 4 files changed, 407 insertions(+) create mode 100644 test_parallel_performance/PIVlab_commandline_parallel_speed.m create mode 100644 test_parallel_performance/Williams_results.jpg create mode 100644 test_parallel_performance/do_parallel_performance_analysis.m create mode 100644 test_parallel_performance/generate_images.m 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 0000000000000000000000000000000000000000..15caf99a2e337e49b96260c5e6a1400d23aa23cc GIT binary patch literal 40541 zcmeFa2Rzk%|3Cf_S&;^1<;aK-+1t^uLu7@LR8~m#JW7a+5Q;eV9@*K6%v4QuIp_WQyw~UTe!k|xbYq5rW3rMmk^mMK0APXt01O%+ zm%pTP?b2mhXV5nySzH{cwJ8ZNn2KQ;@buA6l%^n+Du+hqx>ak<` zfGYq#E-oJKVSGG1JOToILLyRPq9aF$C{COtAw5G0p*cfIO-;+d&O}ScMo&%6%zvKk zA}0?I4}?ijRDer_otuYij|mn50RhnwB64D4a;~$~XSx3H59R}K93T7TAwO&^I^fW8 zEbQZ0m`VTw4)8G6KYoC(KUjybaSr3+;S&%Z0Y8v?3^;^^jeQ6Q`|x2L9Pral;QIj1 z@xv$1a$dqcsiccXXHLrX#6Jq3{&HqHnR3fA1NYs>PYDQ5ohGNCWIV^je4d4emycgS zQ0U542}vnw8QB{us%q*QH#PO{-G88OU}$7vX=QC=YiIBJ%+1~7xu;h^;H#kEkk@ZQ zqhn&@;u8{+lC!dNa`W=v7ZiT1sQgq_T~k}v+ScCD+11_CJ2E;pJ~25pJu|zqy0*Tt zxwXBsyEiT@0Q<|Zz~5g6_S3kIgX21cgM*EOw>K`VLpITLe_#Kdga6314N4XVF1^gupJt|T-yY+%Ptg6d>;dh+xK*fYQ>AK#N631eT}Gc z;VkVKh3{fdV}Nbm5Defy{o}pU!d*KG^_D34NYrNhiVd)-{vTm3Sj!7PUY*=GQ!UY2 zJTn|@rJG!v!)8MH64tUisaoBP=KE>hX2%QdWzu_fPTy(nZ2+{FujE{PyX(j1Y-RJ` ziT(@s-&?<%E7YNBRd0E2u{J)%Pk11d8YWxQ%8ZSIq~Cprk{U7oX@jOqB*fY;3@lq; z;7z?q)&$VYUTR8OBbiXQ{(o#o|E&+YV7nFX@&>5- zQ4YoYw5)kwm+y-(izG7pQOJq?DAeZkzdVFK|NEHVOHO160JKkNtYZMtxf2+mQY~-C zw~xEtr+~=Xuexkxr;}`jCs8UiHuaW}MdAYK2_4&D299%q(b2wgd}*}U7=Uq3&`~u> zxE9-R)0jzmaByY{2SsFvGLp_2P}DSMT=yZOUreG#WADN{b73QA+R&eMk34!a$H31e zaTWtSIg7VUROG^s&ws=1nG6G-+1AkE5HC}{fIC?w+seY z4;amxJE@HUURJxq<{2=+u8g^sC7>ih`dcq73Jl=)3#y-s4`#sClk^caj&`8V&k z*AE;$v={{S3H=~8u^QZ$cP^%Rr<}>@#?Kf3;@Fb_sJNgAVXm&41jv2^6IBP zT`%ic;Mt&)8elT$3d8`Nsk5qX%xm&JWaM=UEfTwR?Ysy%>S(x-ZdIiS_M$PvfYI&@ zt7(Kp=J3l$&{Ip>YVqMC2rV`dcv(n?2lg~;Uw~y3wx8q!U?V} z33?O*aQ_d>YmV9}8^QoR5;98PE!@{$UGE<#*FXKnoTaKHyB&pR|3*f_LZ*lK+<-r? zYH|C8hbJvp-2_IY*qoofTeL!?%*f}39Jz-9Xt;)k#o}YLmm8y?D7EfNDyrf8<}N7UIwc!}5;i_OHe zr_7e8c01dZq;T#E_$%Q_C_Pb3P6$chFnsQ~y4GWVR_<8qM0|nzC@q#a0N??Hn#>u} zlluw|rwA1WpyLPUWS+_d14IXd&r5mPD!nZXU7k0zd*0aef7Amwd7E~KZ1vF}6)`|; z2wA>H=oHdj1blXMnJ_@BoL67uAus0NeBn1C9hne5$WJm#@a=8IEvg-_fZx3(J!EAW z11xf^@iIz6$3`$fzs{tQEWjM?^_wqWXc<2O10a}_SOn0>p(Svz(Gi~$<6gU6=aLaP zxk*8CH^gIkLR2)-Q&K0km_ec8ho$_zOssn`(ny223q<^q}(?{(=EgI;7 z!o@}1OnWu7G~)YgpYG->JIANC`W~`^i?M4D%8lqu^543o>6H&xrgG4@K)*Wdsiz|@ zgmN-6V$P7|R8sTG5)3o>h5W;1a)pgcf)j#rH7Q$H?i{#l%a6WlMgqz|zCQ3dzhKh* zTcUFR@BjHClS~<*oYyusAdIYDjGXWr(MAjyIn!Eu4LF9G`ZVH=I-m}l-{owdF4KqD zI|YLXfq16*RTp=Q4G)LoC6jw?3_7K+^05@_1kT7`*AmsyzO$upC!gu8^ateWvN28j z2Oh@vK7WqB>HqGGhXU5Hb(;OuZk(gS`=pyN7L<-MD(AeRjPU^4jf(hOn*xO=zs#WI zMA>eKy%%e1f-GrF{o#8SbqWfrx)P+S0ly}$ez7E-KQG#EO1OP?ca7nzBc zcwg}54hDFWYM}DUE;`+(mxVduX0wPqRni&NRqh$ZdJJG?Z@-01iR0;L2~;~H!Owx} zn3)?|jrn?w$4-o&zyO6Mq7G;kPzn4(0JG-{V2c@ww=qBv#XbQDd70`UQrf!8CM3UJ z#g9gZg`j8GHWYu6IX?ytoP)E&W6cQsZ_4a@3;EPyYbI5B>pCgE+0Wf zp;a+}S_xttq4MT0nT`LU8#gt5%MN=Qv=cnQUQ8N>Slm?0Q1G^!6{Z^@knd%gj18zJ zdLVO42~SPWg*BpNg*3Z<(}bPRM#L`cv0^jq#`SaXq&eSn#TKNZs?yct&?J;`WLBp8 zM-KKuMxBsgCbsJSG32g#_yVUu(+J|_(k47XFD5t$@?tQ=)Lp{Q<4Emhxa<{%XXDjx z(fTEp_%A4NBnV?7m0Xzjw^QOik&;jHEY3qMzSQZan{0OC(mB3w5}CYyDQx0)Y00La zQQ+6zmLAC2e&aC)7`Po#f&o~#4k*Biw7+aiJs|^D=l%8S^2^N{iKisLf9HcPBx3y; zmdRULI&r0QRK|vrW8)Dc(K@i#^pI8;?&81t)U;oc$ts!DIBm1l<+{^x#sYnX21KIo zK-@RO_{pSQC4ulou+iM-#M)1T?A6>H?u#Zm4Q=T)x;?&aJH9y`og4CkPmNA}>WUxb zgMWTxoEUqbC{E)L2@cK1RB`SNn+o~z%>k#+(;5*bEbL~V8*;i**=v+IUQKP^8{3%C zyz>+bk9I3HEaPam_INV6Rwwt99m~|F`7}K_;nak+>o9x)5MU{imS(<)LK;}GuV4V$ zO-ub{*`{m|J|#sLB#J09C{7IZW?h_Bxr3I>Td-mkZWi%u2}2~f?oiG|`Brqbp6JSM zdD9-HOqrt@?Zui}`TdFg4Zj6_2f&$uXs_RV9r_NxsT?5i?{LO>iVo83cdFC@3P8-A zQz^G~`nn}1wWUvsl0Md&W~D(?+6B(gtj z3Ir$)U4we=!EI$`y{DYi(uuCoi60U=T0Ssx7~fiJYPvt|CUPbJ0`EOgwA>*%k zn6udysOu;xt9*l+zo|UzDHi0vRFHm#;=kMeA3T1BX{qu*?YaZxfODg-;&IqHPX--H zid3zPBU-h=c|lQTSJp^Y(so%xOHa7LMkO%-e9%xKX9meO@xpjlUAzCPFlq%&m#(=T z-+ApqrP&dC)CiilNM&)oNVxyauyXZwHStatQ~L@jY+&}5b#Tg3e_!7X_zL}+#B+W3 z5&jA~>C$8iWp!St`-3N+3!U1|pCR)lnLO{R-jP#9;AA#Qi>3p&fZtKjRB50P{5sTB z%vc4YegjRKF_A6g#mdrFdFr$VF;f6K;1N_TC=ve-;mb?Y5;n>i3RxwB=X3>RkG86v z@?uVtriZR9!x#D2$QZvV6u*fj_jau}tQyo!-f6&P3q0{>J{dHaV# zi}P)}{fcJu>w4GU0A@0OQXqc0bMK@u1y&TK>dfilq4Q?)F-Mq);_K;YzeD4sm#)(I%#ozJJ`WyHJ3_F?M%A~)ITju zUb%4#b!}D)kyeiM(Fny)r|;o0rrW*Ppm%HgcIS1?o?g6pu~iJ zt}b1Up9@0c)_$`Qum+t9MR`A=4QT+*lwG8!hB+h{!TKB}AcMIB7~tU=G6@41YHh!Y zg~xRba-`^t1i-ixq6_>M25v()3T@C=%wrjA?f-@(OLof7p=envFpkjlCIR%)Ki2u~X+>UvySTrM;fj+l&q>B+LEI zoUR!g$e2HJ&S31~YHu{6<*1KDu};5EbV?d34Au9Zr1jDbw&x(+>nQc$&C#cLxOnA- zR#dx^gvQgaN*exw4MWLqdmFZOdM7j~(P*>q@u=hwXo<)eo4MJZLZ~zqmc_#PDm?l{ zr5jEFKpkLw*xDBejt0X#-H{+_n(Tb2FTu zeZU!~eqx2}{t0a8kVw`T4~y!ny09D9$RQf7a9hnm%^bB@P#> zH9wE?raS7Q@=E;y6`mYPs+F1b)&pKv*?5FCxT_1j_eZd=ZsLz2)?j5AAj>Nt#&4+s zjRDYgEWfWQRxX<%pASO7lh-64?WNh5s)&mLYJ5N~6M@?Yn~qB>{9VDYCmeneZ-{RBxNGdR{^9!GpXf@ldVn zV%`bMfe+)sBXS_I!KxgwYIBncSW(_E(c(M!($|= zcIe#J7L9tA@_Jk6`>D6?Z~E$?kEe8FfDR*2zLf_1&4*ia6u{ow+clshEZ0>)!@3opE24DA9+8b%R5I9-)M~i@Ws%%Szf0`uK`>s%?_c#;UcALo3A}Ba&y<^ z#!pn@3~uMvdb85NHom`h>~cFphVk>6QMz8zpc z8y$Z%m;J{hlYJY;_g5T!wAe{h75Z*$;(=}Hmq*I?Gr2G40%E@DbH6*4-;B0F-=^qy zZN&fA?KI^1q3H>ze7>-F!_5%fBen_ti)NIN4Ies}Lsg*VLB*miG4@A)Omr7oDqJ13 zGV1K|&E`#GpZv)`^OMuLFWTi|@YpAoMp))L?DZ7-a@H;u1`tEHFa7AZ0YAt2lketd zhjbzH|9G?1TgIOXAs$oGysRUnjqaCDTwQ&Csr+#w!4Y)Jhe(_sNAk00DXw^h?}ctm zd5X5tu@_W6(F8ZR>^>D>9eSz(7bi&2uzexjzW&BaPHQ`~?#xF+t@9IP6Rz&%%6Qks zBk_N~>v1iPT2b3iN9CNjkn&u#FRNOLA6|cJc76C1=;Jt6vz3+%RQ`V7@}?t+6n~tt zw8S|vy4IYc63Hpa*3YNj6|drN$Q&V4WqU;Y9%tCzYlih?l@4eTc14%xLwC<7Vw-|U zf#??p`Bt35{&XFbBOuiHe{cT`pBz# zdIz6bi&eQW|78C9pAqBtkN4jvD6u?^L04pB%ANy+14U(}eZQ;1#M+@CPVv}+y7iqR zxQl{$XOXbZ4vm!#0=EY?-Iv^#me>_ewQ6-S&BM*uVqGKyP38+?lq(RlDARTe1|Y1~ z!2nTypobt3Y23k_hAz`7MEHAbS{-d#(htG_A2Yi!z_cEyj>lpEO6Z8j9O&i(^`Nm< zn2&uNsgrY=)fO$fF{%eOzgo7fCIniBiD297@G&VQ>lql@9g4J9T**+0_RHv6EFFjR zvweGbPztFHUML+*o+QrmOn6+7M50vS5MNQ3|$#T~*e@J+t}vifAFl(uLj! zR>zX0Kf2#a*?Jw}i~$anTrOz1@$vZJY@1QAP+m;zK!4S{0nrzcCzGd4hi5-r#Br&> zQRkHnNMbiMwNj3wR2uiKCL8ISw$eGaP2mp5g%hb22KQ;eI^sSGRUt%6cW6Ev!M)#w zCipBgkG9%b>g4*mn7v9t8zqwD+O*kH$)T%T=h^NRv<=DL-Bed=6zP8blsw>r`#|9P zW0x9)hRTIB4SFNOZMP_&bjTRbg@w>OY;eJ;y4!l|ypyHxM3M=Ym*e5sQ#LM3dGcnE+Pcki z@Z!_UV`);AWxs>=-Y~s);;4W(P80c#=b6Kn)@VXsZkV}PZnE~QdhpOaA=l!&26r2+ ztH=&zh)kKvOtr@6FW(N-vxw2wAV`gE)#Y2UBf<+iJ{xfF%m;W_=0KwXc66#D1+2zr zlmvm>5ALXV(Dk+w0)k)KJ@DIv_A&zcK$O>YW8DxGfED|*yDGt8TR{OL`1(ZuTr8L{}ZZGux8jxr3% zlrmQJxq&u*4_sQ6CLe%GMt8|xI$7A}LS*ao=0&GdW|;sCWDotPG-%vrhZ9FBs<4kpA6 zMYHmy=VOlgywG~eMS}N|WJHrs6W1qRxWwIG*P!uY!{IDR4SQ-rpd~tH74)7*C}Dsh zIzKwd&IAOt%6v$c1-egvd)tNK-~3_V0CTEeyt4TW)ebW+duHzySmS}VlfnRwUO<4C z-@$)Gp2MbI!M185378X+U+I{GcmJIq1-_YXv@enLLbjAZ$NlkHBgOqW{kP9A*9F+P zT?cHQ3UsT)n5LM%OKrJPP5DlO5Ssn)JvL;ub+(|$56k$Zc0J3SgKGaf?J^)D+?XZ% z?S*=(u;w_07{|4V9qn~W*2c0C9b=vKGbkOT0S0h$N`p1qY`|7hMN7) z2o33Pl87#zLz8XhVt`9S9T=d>Z(F^$3xx+$2M?65ABL!c2Ck-wE$$YB4Dbj!U>5U3 zJ2Ne`Q&|aK$0H1|w&V19UGvilhu1ospqp_p$vpeBIknO`qJT_VU3GfV$4)FEZb*Ju zRDZvmO`=oMJ@I6&q4P~8dPw!wRI2h_-{j7h;o?U$P7SBu(nH9eiQJ^nYfH0idlg4j zd-JVT6uDyXiOE}y@LMa45iTZl!Bk->Mb+m`3=L^^?z3=pQviSG~^s`q%5sBS3W<{!`se#D#!m=QCIDOjgXZ01>Rfo zmNO4|&+siudoN4C0PiMp+Edp~2-MUF=tcxxcf8yx-wk=t4$W?8>cLGKhef}!dDvGb zShW;l%{Nt|x?(-G(&T!ge{}A1$$k1P)SD+vedda-1^Nbx^EFW~E5*a669Z6fyL?;h&U@0T+?OsDJ?Es5C;-%gANH1a(V zi}eo_Kfe8w3;Lh;ciz*Qc9!y5;%&+?MaZXGqZK}zbZbmm+chq}Ubv|)M%HnsotGl# zjJYMkNe0>XJff@N{8@5CdSQ10z^aK}tfNy&oYd0fT;JuEG@hc15+ymhr_F1{ViAD> z&bjt`r`b=jKvs?0(Tl=VLr1k`v$@X2_8j}TTE>@;?Iy0VF)upjbYETkgBal=#pR=7 z&mK7yHe*l6FVXn4K`mQ~?fI{v+V&iqUnt3y=>hI7`90_77mBJn#}?28HFT=4)7chzWEkribl`PWF^L{o+T7z?}Y~>*p=1-#;#`%u&%XbVWOW2bB|{*#2k!$G39j)LI`}HSozpm0TNhU#5sS z^gLE|I)A7=91Ky|8r}8OF-Uikik&Dk*ut{L+S@T_-dm`T^LNPXMWc0O_b# ziAP(a?s~~Iu^hH)zm7|9d*zY|&6bP^mrNKMnh?yf=vP?b|Chtta2Yd@enJ^2*N0{6CLIz?Zg&DCEi7o+WDwDDAT;ow zQr%MoLeM^iX@yf6ZZXfZQmDFH!*8~2-}1FVA7_S6%_<_sL2!!y3bMArgj!{djxGR= z?eAc*?E_4Zcb2{v$>B z?2a<_Pi3_3XbfL7IjMbd9{4~V=uJ3h{Hk4yYFk~UBRyBHh@f|K+Lb5C$(&{t7mslG z5w4K>Ey{~AM`ev_sw`d&@#J&fUzSLVbw1G?v%=HKe+qZfa&=5mXrV{(^|aR6NB8u$ zw6NY+af9cFXZi-KJhamJOlf;u`#ev5IM>vII+kOHmrQ%kV={o+g!nDCYlA&QDM2%1 z%K}ffA@kGPUO9EGq^r%D3;|h-gGXhS7Ftugx78)ua*ERUsK^|Q8>*%i7Fe}Uku&zx zK47p9Z4boz)CaE?Ba*Xtk-@d`%C~;=0S?^c6I6C0JPMjk$Kv|Fo9&Vm@zl-}^NlsA z=K6)xA1T>gt-I7NUNa@SEn9DH!LUGWS96uU+6Y1T%=UCeN;2pseZjPLxVrZl)F80t zB&QqmVU-}{M|6}xK{bikcr z8!6n9YvC_tm87R8aGtcQM{Sp1-ae_LN{)T)Jpk)F7iWfVa3N-JhT{%RWWdWWPkGr>j9gsXJKINZYVJDFX^)ffe6^)D-!+I z{H4_IgC5+Im&rQoHsNeyibv>@IoIvbxm%b1#Kn=UpX(X56aFy5D%I&{uR6AkpX3a@ z(BO;cNl+3e{zApruA9C=#D_}&{9h=av6e#p#ZW6!VBibI|Br@-?brV}3oHWj&1WNJ z9hr!Pc9pO%h~yit_!A8Oiaz0+Z+d`zN&@@T{(^bUNIE*Hywe%EjNst79vxZLUA23;iqjewqMz0BB0(8oy#F+g)}tYs6@ zbqoWHg8JRFD<8|&FBqQT219YS!PKkJ0nkAHpTp2l{iE1vHKJ0g)=Sur=st=;R|azn z4t93R;82CG+0?ZVwzuYoa{~FH^U;C3s>y%O+7w+-;NIA2obu)r^Vus50h}W$L;*0u zu>~>4nXRk=HX*SbwIk5HD*1HcO0B2vJ{G}vf_T1Ikw-hd65?j9?K|I?C3assm-}9R z<^7_l=yO3FYr6CUU(SDia4cx){^#-uest~q953|CoD0jJ@)1a}UW_hK>KwCGRBnw7 ztPdSMk(u|1|EM@q^GmGFbMMJk*fU9d-BX~hyhEav@OL@KwHCdbCKpG=t#@Dx&Ot-A zg>%Z8%f3lgG4Y7)$Oq4{0z{?H)KJjc&=sUWaIQG%9d5rZwj%LS_Q|t%vv=qbPjZGG z9=HfWaA$1nrQ8-w%vzZjd*UJtMiBlv0o<=>Dy#*@yatnn#3Fdo*s0Ds*x1t#p#AMQ z-mp%1?u20rG~&G@vp;aR3l5QLN;x@<)Z1y^mue3?M-jm^+Bv*jbLD8bU#?lI>~qr7 zA8y64%T$98b6Y}G`xz1WOdT@ zJVqz!m0Sc60!L=ne+AB3Gi|BI4~Ca8cf$v6huK@N|3mtHrDN{Gtn5u@`darS;6Cl+ z31%1MCz#3Z3DSL6>-@y1)b)h&_k^fb`s1*FjE+O^*!4Cm-zo3 zUB*aPncFi| z$e6#u_c=+?dt{0!aa}|`vf~21q6w>y-xLT>E|i8C&@#4nFzn5WJ=Jq_K?N87c=hZ<)JWUhZ7t(iEsYPfIn8K# zt&c)dq~)Z}Kv&XslUhMD@=UUknt|pdK7k*$u*{gXuv;UHa6H|FmsqQ^Q2*JD%Fx$O z>&i1SG8XiY*9#ZvlD>m5Q}pAB9G@lztZ7ata#Wl$o-i`kgViSU%{^JmAA;8Rk>zSH zJ;^;1F87fyaNyR9TGdI91@SuynRPW8+itjpECQGW$=mx~zeU5|i9>U*so_#(uq3 z{6Tn*3l)f%%nSANCnU#ERWlJoL3?oRD@oYJZ4Ta>+t_XCC^VkW<6h3L=s6;?OoC5oOc?=6U~dV>%W*{dV_EFEMH`K6#T< zeA4nh(HZl);bn_9ymaG6nhUIQEeknkm_B65IAln+_0RY@>RGyx)ZKd-pf2hzj^p3u zI_*r~Lbn3(Fw610$C=hHf#agDXQfsp`AA$~$!pRQh5d(S6)JQo?6ZM1!W&9vZgfpf- z%GP=C<0xnarY+8H@HX@BFFGd5Rh4@X8S_iG<4;NN2aoypIZ^xc-;o>j%ZR*xPq2LX zO4`O|*l^nIHF>;+M}6G%M-zp?=hX6k8R2ja1|YJ1z(L-BQ{sG>Z0Yh8+`@{n?W7>L zcU@BiK`%iSUaq}bNs5+%CfQmmg6=(=&+_FUD?^teK-IA2yzNg9n{H5I)@l?#>~fg|Ax z#5xsuHrAa!_i5UW=iVp0bB~eYq9~z%^CKa`&Tnv8r^^17{gr8()${%)OP{@8v?o1P z5!gLVmJqw`C)*iE6|_4F;;h7GibapKRH=IOn3B+CTERT{ zS~K6bOc~3R*Z6pzhwF8~G$)his8?@{uRhQ)mtO)jX#n$j5LIRl!pAnE8)+xM{qfOQyFp zV;*}kr@MOZZ>I0B|CcsaRNBm`-AzxsxLf+Qh0m*`XA`nX9zBw3EuP+nSLXRzW{;8J z2ws_RN;iD1SWWc?tEa3l-Z*Earpyd`j{e%fW?5F=K$(Ntir$&c#p+U4cuvTO8Ipcw zUX13Hv-r;aRb;>3_LChjOK+@aEAnV|MIn~~C0Sl*$O-M&6+D5pTXEKR8>?qOE5fRB zSwprh8EA0YY>qrR#efh1x3>oE8dY5r65CWX^VT;#6O0OzODI_$?U5(q$5S>5kUt)T zw|k-Lan!2_?tYl;kP@@_ndZFy6QnL$T@*75yPaa7d+mqPCew!m=}t0k5<0N?kV)EDl~yq?W(0c{|ws2yzMN!?;-deIfom^(H^@*ftSczKPQIWY9-z( z4!`>zT8+OWvvnID%o^6C7KaqPWeuCOc5vi{jY$P<=_Cs|LbpN^5@3zqU$4c@u#dNf z2k~}Io<;~2J^`)l9INRbQrk~rXC>|yoU5~_VidgW6oxccC(@>xfEOY^&)#8?uxQ&F z8-V$D=wg845>a8a>Y`KCKG8|NU1`nI;RsqYPKzidMcMn~D%+uI^R)<|>He5ZA`~QjSCWYuFNpq5=j8 z|CR>MjP-v1kL3Hr9{T&b9FkwJd$`U6y}G`>wA1qS+6KbEkzec(5GMXnh8+&wOM^ou zJ`Fi1&DWYtTQkb)&M3~QCH|RMcOsqBayEF+a%ED@4u^lXo}UrCjNq&tK4gFI5y|bK z4>RbC+ZHxN>V@u}q(V>ft z#V?80nfG%l{Nz1v?dOiXr$88;S@62hw@A*H-@0Z@tR1UkMkqJFApS~Xx3FkVUmo|_ z%a1-O;v5phP4Cq4u<{=zCLB^i@C0X^8xJ1566|&LZhGui#d+|4Gn%gsRIs0b$xq2_ z2dX4l)0`M_D{pwK2t@zA$WAjd7H-}QK7Z59iQxhbnzECxW*c7>M9w#Hqw81nU3=h) zuLgUL$}+QYE{02uXi-8dZ;8!6w_@(Cb9O7NQtV1Z!h{xBsWg^L1Jm4%Q|@H@SFRv4 zcD9ff+GXlVa@$;O+1-n~1C+aWbu-p3KpL!gqeBHpZ|eFp__(XKJi(im5G(FfB-#%7 zxW==~*rwe#b&hw|InfrX)L>w!b9~rPtSPlmG&XX$#!J>OaTcnXINIu%tImoPTh(`Y z6(Wk)Se%dWo%M3>TpZjy9MNez=~*c95I-p;O?cK9H`H1eOdO*FU2R{>L8%<17##qE z{iMV1ufLK?`9IXQ`Ie^mt!B-EIQuTtA%)5@r)r>T~Hx8$$GnzLbZjVLSlIAX}txnsdJ~T?XtJg4@ z16!8rHd51kOJf!E*`m5XYTQXmX}cqNtsVmaLxv%^TNU}OtiD&T^=VC0MbT6RG*yQN z2MT1bkUb0YK<23SVgPMRsM`_J;RQ{;wP}(=$?Ok#3X1deBO&`Wjk0+b%oI9l1ldeu zLdGA)KdNs_jxbgz9Ajj(Y$UjR-7!lI++#NwzQ)%ifN_)wRDWG?&B2dDRw=X@m+x^6oKGdmU^O7Ljx!NdW*ls&Z;;bN;=w?4d+qWVxJ09MY}N-kUWwyg)9 z{C>n*jS`5YspA;HH=k_!uwU=VakGtw4)&mp^ee{+TNd?V%XNKA3ziv8c^SKu2pbGQ zK;pArwvB2fcQZ!i+JBHxMsl7|@&iVeKGajU|Xy5MgbhNd1LAJlAnk#6A zj@6mEb(1V+Sp~5H>3p|_W@KEozPKP_H-d=dwr=J`Q^e43P(ad5HD~$MV^F6S3Bc6r zw@=D~?hb4-i5HE}4BpTablXIeH$r7A3b*eYJ?^yWI?l(1{kEk=%;TPH4sb6V8Sth+T_5#J zTSY!HULi?54;LK@)>vmEM=z0(z` za{827{&VHYUW!+n{@5^~+aGhB3H_!L(YYS*mGc~dnXZOoLnH%tc>Bbww;y=>wcm5j z%{5Lo+=@M-HD~9!KphFM%Q8iR1Idzkr}28r%}Z~Y+BwM6YPqz**Myg1L{#ea6ATc{ zJo!n=Hp$*_$KxEXgRIhqI?{LeTSCla{jk~bs|!Rz!y4^}IC);ob*sRapV`9Et{7m4 zMpGKsKSl5USG7vT5Pnvh9|IIapKvR^)CvO^d6lrHj8AwS z%kVX@!XTr5l$|GP136P>upeGN3>{HpLLCPChXodqV6j26bVH|EnNT~Bv41o!+iMK_ z()b~V34MUc*7b}Aq%szgF3JfeXHqGO3O%6feSkXNkyg{`?_1LyZ{Kb+{-V~};Ye)4acqCK`pqtrX zJ>ZOdj8ycnKjbvM2tmy#u7A4e_o0;98&)I0xnT%;4m3GlIMJ0Z=V4g|BD?||C7Z81 z?v;NVc!88o&RBuUVvlE+m~B1vL%r_YwY_o48(FBD32MarFrl2#b3qh!hQV^^g+jj> zsdML-7)mcKTx2znK{oO!dZ`UNFUba#4ZbFsHsaX6son;w&hvyT4ua5{k!<>SZ(pSa z_3?AiLKl@B2kI+>3B#^)2bieLG_x*+~;i{^Q{~85NWBk$v_0hD9+eaDd#w>={{r^KSMIw zFDapT9IcCr?5a>;9k**j(%t5QqpVAE?48nY>|%N8Y^xVm4y~-DPl-g_DnLqSNJRT3 zw$9nC>7o*DIr8Z+*o``}X-~@E3hO$x4VDK@{^-cjR=wMBd8TZD5E_-f=oR;*A|rdp zbJ8|p69aIpaw3^n6^93oKKNkBvrYY^`6Oj%_b&^K%*fW&jF{wGgY}^}f3Juv`9EL3 z>HF#k`}BXSB>Z|2;Qyckzz4B4Uk=4T?myYr>i4gI>%qEzll|-8P;2*_?R38`=HB)H zmx@&Xp85g@bNQ{x$iK5h9QBeHTbk98p`)6EeOD|@cLUA53@T#og63h@@^S2KwcUyn zv{!w|RILG0;B?!J#_rFBe{ts(7c=Bbgzc3wp@k)U>kBh`2E2V*vJ$Kb+SisJw`b_v zr>XYib@m*+`}D)FZ1!Be{}RhXCplD3DuRIu!!n_9Q7 zR4-0>-m;x3Cr?qi{(rqx&pxO){$0Glw}M6orr>?q^ZoVDIPt%4zeoS~5c-XZ_Wupf z`qe63|HZin$te5Seee3;#FSsh@BbS_44LAMj)#^Gd`0c-%J}Nund5MAp}_g`v3*%5 zjacS0555s3Xi9#?Ldy4Mp2lEda9tj$zKf_+QojzXWF~rzQqC>c_d;h_{}wK_s<`zKm|JfI%niO~@vu#ux;5ET!lN z$4OAZubsKPS0V#M7eVs`A*^Q`K7mxma|1J_j7oN znGP{{>?&;#qCy5;5t1G@TihMY|G0tAX#9u>ZxhP=73jtWZPecIMO*y@o%Nf38!@wF z#ii?zUQ$Q&q1}lqp#msV}4NxN++b}jT*yiWz zpD==FG{nu7`pi5Cr-Z(GZM9G4x%zs!!YOW~O47Q&QB3kb|LHt?_7Z~&9z$sP%4u&N zIB}Ik`{4aXuDllNg7v1zEohqEPS_{-iTmftPf?1CY~VY2&YeYx1$D?BIiGIlkh$0> z+x@Mg+&Xg2)gZBJ^gNOZR&Rc&~b z0)t%OgL=(JpNkf*G+U)NRX7v2I`Vlff0`Q9n=E*Bzib#f08y+Aot~5ObYrad(BUsb zV)MF%>o+`yG3Fyr($rVjGsc<_38_Y2eaJrTO5<8k%lpVs$mnf4+<+v96-{>cpr$;l zPotu>HOnI{y2Q4JDiaUW$~Qe>A0zUtQ?^W^AXmvmNekeQ3)Js&;>cdiq@}4lKYW-R zz?wd95G%Y7j=sPC`D)uLHYZ+Rhq}*4$ypHo(%aY)E;qmIoCHqU0x~Z?(X4ns0*zaS5j1jIHN^#PpX!79 zN$S>$&esS`Xh5o5h<{3A?m9cRpkP8bBw3c{x^S{Q(Mz;Bv* zv6-z1mg$f)Hdv|5Y?!oIpw)8gS9;Lk*pMwO4b3HV6kG>PyJw3DQ;Mk00vmfMw31;+nIu( zHenxQqh(n{z{go+3^ar`Z7aPgn|ce|>H(|i?bVk0l?HqDvOW*3H9`^dprE1%L)(D~ zzfSPwy#iHh{Kz#wGy*yz`K=FoWw(Cg1|<96U!z(Zwm^V*4}M=1#Rg5sxqqth-`=}F zc!TfDV^NfX5gkQ-iySNP8w<(@bNv@iU4WN&PUO8dgLFLWAvPJ*kCivse!1p&Gdt+`s${7ABFC?xOA4*@F zB<(F#_#7YdCK4a{r@u8rkcwo7rMC&?dZBEcKedFal2YH*)bUD!vCW}Yodt>vlMc|m zCGB;MYX>$cST}4#Hm%YraX0SOi3D2zRu{oCk6x?BH^bp#R5jlGTy}xAR908$vt67% z++mIXxn#Gh&y|Wov6P0*FcWUe!Rxr7MDgI&0i7a0e23f^(`nnHlgw4Jt?B~AQX_mX z9;ZGuWoH|o0^po@!*(cRdi1d&s%auE$qsbO2~_)uG3YYh4oq{yh&5SYwpLSNuIWaA z2H}N4)}S2&JL#3J7d*5n>C>h6?yRv&-ep4LaFxZG&J2-R9?3d2uYHq0C4K|h#k8J2 z?(IzzpD>D1yLBt4I81zlttj(!L!9%J?p@%L?&8a}Moat??EDUYC-uaE( z#3_aqjzNPU>M_ufcEmbKYNyAFVVP>iuzd;HM5d!shh@(*%FSnbrfq=-EI)Dso+Dmxao>$hZbsuk~DC zap~%a$0UoO!i<*_s~fYC~_hE90XLsg~zO8xNZcHB~w+Xig!xKdJ@P(%*B6tPZB{P=B`teVHPJ zS3%oGELGV>Vq-aF?(vz88{KE6x9en_C7{u!-R)tRp0Pjozow*^*iYi{IGK-1QQm z(}?n+_rWDQE>0*4knVcDPm&dR)y3w((O#$_U3I%PW-`QP7lbL=33tBw)BjERQ{O|Q zJsPwU(*{g#RlNLU!yA(q+Z}q$tqY9sSt+?`H&Y{O?4%R~a1|YU&n8|gR(!N9Ay;qi z$M3Rop#mfSoeK5xrzG1@n=`&c;d}PruOyY?R$siZv*n&O_$z4w9Ws5rumi?@;>YLv zj-K(6RzvRk z#nO}gy~6U*#*Zb=iK)IKgKmNGTW*)89A%eRDq?0YazHKR;gyaXv+L^0)$%h+=j~&H zhR-M1n^d<%=-yw^uY)mZv`)FJQfYsb!s{}=bfKa{ii@K7Wrz65_neNq7+UYUN6%`?Iq2%$|7J~!=O)c5GN)2uiE9UDRM*>eP z8lLA6rRa)bOQaFZkKrJSJ?1pi_hr~t6B)~0I<1g~W~;WhQ-%+oZH+y>gfzcKL%;1h zg7lG7;D4EMxcAHllC2nYj*~JsnzRiGxW*?TZ1G|&)@u&g#~ zEFNqdU5r>gU)shVbzeE5_P%+9sP(Q|F_hGi^Zd>t`&kU&d zWT0;xzTVv@0kP0_rcCI>RsSoJv$GCKH3@1DT&>|t`PZzG)4j0L0nWa)y^OT4WX#ug zdqLbkvDnK<`^t^)OZ&ZO>|e4v<_3oRgqL+BF3j2`?SJw+L;Av(jn%Th9dFuCd0Lho z?b|7ItC7@w zWC!Q5lfqd#Z~cGmU3XAZ`?e0!L=II!dXO%?_bLhkB1H`#QUnzO1Pl;*Qxs{T9#9e> zNQV%5=%CV(4xvaBL4kmjNDD3J<-D1D-}Rg~ch20IJM-qf_uu|wW$m4{X05%}{(b9P z$>T)ik~A`xroKAKVIlF?T-oVY6G-0BSNI2msFy7!%g7fvO+A<}Cfz$Z?L*HA-AJez zx%z04+c`1bIp~kkZE53)$+u@Fx7n?2EvN$61iJH1UIBH*uFgo;!*#2W=~qQRDyAA= z=MjTVM<4_0u!a{bQDJ4PxT>$@GiSH4zgzNwWol|m)~9)nV)PB0KU-BG>W_c__oJG> znc|;D<~a`BHgeUoVe@Qs(CBZ!F#WGu=C1`a&A)vA|0Xm1P1!{7tW7!Qv;jtsR9QZ( zlzjY-u=ns^W}5jKXZh!;uDB&T_CuSx8q(0j?TkP$bhQh507C+seFfclyoHpwJ_VGR zQYEubg8>!6k)7(&r-!P1D;oo3{}a{ z(IFCZ7-s%&${8+CDJSsvVCJjhcvPdTO>S=sQ+hwAYL!%dz=_Iv^6j=6&#O(T=DY@9eAbfVZ zXD&Jviyfyd7~`U3=GydquR`LJ*T!EpMN%)F=9Dp;uV^`a;liwa~bn%dBbDOyiP9l)q} z3R}5VQI#R)RVxzjYRQ-)Qvq`G^w}wC;~i^R0Mb}txm+r8l|)qf3yLp=LUm#@9eTC% zr#dW6R-nM0>JU@+#QK1zOb+u#z}-O;rf#u-%UqJthIXlu$82n$LjUlv7A6G&vCtT- z9jmT^?6C%_2N|o?#(R3VZKY$K_;P3BWK4FE*qGBxi*yfX+5REHbVh#+iU+c`NXd0L zv>1{Qj^F;Q{O8a0WR+D)QP`S)B(*XYR`3Q6gbAuD#GnEp^!HR2)LBG7%9AF9OWEnw35ek_Q}omD zIWnFtl}xYCN{oD*45{F}>kD-EIhO4Lg7QgR=oJJL2IRmd`U9N({1{uy)wQFayL z{>!#vJloSvv2FY})KCM&XX(LIO_oCYi27wxroIMO}D0tLY<8cQkC9LjUW z;?9s$H9LJnIor)4dKF7;LV3;gohbB|4IlQb64$kXhCo}Rq~??K6GQy?3<7?Et8cfz zdiu4bf#A*Zn}@-t%^WkL9Yz(Nww46D{$Z$~r>tbg@X6*G?W#T8%Zs6)$#%9gS)XGj zI+DV$vMjuch$ck_8Qu!-6)_toX^Z_SA{sZWvRt-xUh5Yo@MW=W$)N7|#kT)RD)A3` zBrf{=6)`c&|DjUn&R(w1FSsR@vMIlopl1yu$o7aBE-r?#1vSQ*UhJ&%R~hdq$G@z2 za`7g&aR?1deh7KoZk%MofE^HtYYd23>u|Jwt>|Li~5tSoW- zW3@*Wo>Wm{QN~@IeonUm7u1AOs2u||I$x8nG)#ugJ9)JAdslUNk`eDL7ihSTPI}aO zJeIKlnBtvB$pY~F{dczQfE)|hE!l}T}6M1cJ@Dd_AqVv^CgOk8%__19UbY2 zq60{GjSKH0q)<621hZzQ8w!6g0mCThV=t>0d{?egfPj*#8jdM&r~-gPZV5Q3I}0-t zT5+o-YWf_Q~_>4cl3rt2;RiZsLV`CDV9zh=1{Xhr4GHjght3F!>8&vEnc@ zn5*m8dPG>B7u1VAWXftu(`VPEh_9BS;D!@ETWW_?6Eqp*SK;pSltTNyjwG$6PK$O% zQGKkh5H4>Vd9UuoBb9@l41rg9BbV0W+?&G?B*DtW5#YUXmoyM*QQ!e$D)5yv3V>khYjnEvxz>(N_koD9J{H{fvU7R zF4tQNJ{oXnAK*-&+t?F`P*<1G0Zl2kUYK{8kjbfq58>pM>4CcHsP+_wV57z{ZqO@B zEJRCB%4dZkE}>DYEh<9MwJb}pb4H}?he%q#D(aTEnu^=&@&I_0A}fl{S>&wZKgsn%zH&%OXA>O* za~SSRNLUf1<8l<;OV>dM75pEC;JdQ>?^pTug}YnxXYdlG!zK3%A|ospdzF=g91o9w z-^o$Ex6QIPgN>g`=YqF(wE!}*fM?c}I1KbxE3wnDFNOmgQP)b#FJeE7Dn05 zy&l2c9EMT}iJRoLP!+W?1#X~2Gmidz$S%c3p9Hfk%)J&sM{r31@OM$S=+84Fs(xna3f#pNnq~m=(E| zC`6E9lktMAo}tf@e*;8PJxkGWYNa*b*@{Cy5bn;)|716iDVh4d{C&`)gzO2=rS7q< zflqVyKPZ(@+>$roJkLur-UWbW^WwY=6I{lzu8L@kUKv%M#dtMcy!N@9MXG9Q96D#7 z>&1A4^y|Fqqj*Fp2BFq}UVyIQWEo@s(`e`~p)R)Ld$VX>4RiTUL`KS`Ox`Lq_0VL< z^ZfW_-Y#zv_wy-3BD@%UpQjL_YK9`L^-AO2nP)m$%rY?^Zp8Wj1bCTbE>^Hzu z(ZT0ubz8HUgs~g~eCV~$L9Z`v_w@SAjp1bBsyXL4$}5$@sOF?VeC*~*l*Ai*&Ltx! zD0UPlU%AJ>UmD~YQ{dV!fG0=Uo1|~2(f3qBpQ`0BQrb6-JSHkpDa-42VcyzuEC~*r z7+86ZsC{l$+n!pet+lOcs$<-$08<2bgB z`m&^O+QFI2YKQ2m8QU8^N&|ib;m!}&t4jf6)F;h}cj7Z$3SB9?YlYf$7(Bl;O(~k# z-n!lNLIvcpn~zM`{0*QgygMGh*W-rw!Z6v0su1A^sNm0lbQXQr>O{q0VRslH^wy48K|o?c{Yt%%b! z1LviY+GD=_7*lwYmsuh&?ne$YMu`esZI*_0J+m7~`4CE(RMX8vxvnsRtMdaj!Afmy zZ1o|YC_gcCfyA6k6um@DZd^|o3*s~vzkbNCm?xc`_TiYs072FuLV|{WhKp zmSL3uO^DX_-SHF#UgckYYTi@mcVd;hU~}J1xdK97R5ZYJn43XaavW)44OaI{bA|Yx za0lD^^d)dSFVW6ttflBC4z1p5tlPl_-K`t5&oZ7bqprCs9&y^~oa-?Z>|4qXVvvp_ zu&#Whh8R_B<;AprBgPoN)bdI92u5su*VM|z`^)$g;gvx(?-o)x*fUkqdVX!{l0V0Q4>t$1pR-4&XR zXW7`vjmkbUcBJegK4=r0-hy~)7q{Vd>A$j6j}Ov%ZolfE zUAGR@ZtQUg!zvf$;)4djd__g~7sPpTviVQKGVr<_+7!4Oq}Rn)Dr;&H3nc48cXlbN zv8Og8*-wqmB+7~*5UDaq%FNh;L;&+cftYjRiH@plNT!)Wviy4%h&C4QcgBu`4#X~gEDY7!Cio70e z*V-f8eMN(@3)3w{^ecg{1}>O!S0!XFRcRVRb4&xN>Yx;bHC^-% zVx!g+;*Y0t*VREN^z0eRvbw(c1<#9Y*+rm5c2pL7#zPv^8Hx>FZ3B@wfW=%7&@z~F zicsSX3N}GS+h5`eN$v}SYutLspH}x{$bY%+yR+7@m5`+J>G|e{nrBpFbH<1A;*aB< zzD8+!8+5HsbL~S^^v8T1$&)Y$NC7J2aDchg^LWt*B=vpMWAh~zTPtfDza*9<9-m@H znP4Pz=1MsTJ72!#si{cj@6VN&+LLi?3XSvf2*2=Uh`E9`o+5r?V{J9D>v*&0)&8Ki z#FVwc+oJVKpi~tfTUs~cd|L+RpdY(7&!}R5T~pf!S#982UXPKBbBzPJY{dRcmLkFC*iNlt8{AeXo9B*1@YLhy|O| zvgDaUAOz7aLyHuStt=^cjj^uCXBQ2HSwEkI{n*!B6a<|pX^$mlS9)kjWo_m)&uZEJ z6az+}*KMD8tG90V`O?awyNe2E>;qk{uC`ay^Xj-LOd*J(4<90knuRLjT2+USOr0;UTuBp#_6KUvaOhpHfTB~|y)|RvEYje7F zy@85||4kZsWUr`bxdcwK)=uVjb+Cq%`sBilu71}g@zNs>jR+pX#+|d+sqX^)-$%}p z3KBb7xGBg?~qhe1nP6;ji7J17&TaY4DuAf;RI#bw6 zjMOIyzEdp zS*YIsa!U0LMv{eWPD-cGy>Eop71Uc3<)XL4>}B=nhYY!LpI{Q9<2)}MOTMnwky#LZ zt;h%fI^}n{>+CEx2hkn zxPJNa#W^-=M(wqpg_F*M#h?*d=nd%hlAXWR1eV_67yN*ngNSu} z@@Mdmn$GhinmG!peVV3!zAF45*%|)b*Gs3g&i*co!N_2w&0^-99ScuYniG8;Ba`Mq zx&MtZC9&Ffuda$@`N6@O-QVuB(=h8hz5~iE+hcA3*gj;GYXX<;Jdan+S{ui(&`Tvy zFvXh=0qIGrdtWJWJkIswEx8ez<0m?ICN?#p6>UhEhmw>UxfT+RkC_v?C4Z+_(u@-0 zZx4*U8KP3qC9dc#oN2taOnp-?$A5<6cW1^Bt|)`gg!cj zYNr;0XKoiz$$q{83-k&07?G6V?9*%BR?~ERf(_wxfw`9v%T40IZ;kd1Qbdw3#9Veg z`ZYUm&dTlvAW&?ic97m~GLOLY3d=pcwf(5Mrllsq=HA57kLfF4o(B*Qhc9q-jd0HCcP}AIXkKQ+t2n1t8`g7v fr_b1Ln|N*{L=q~@$+(e!x7+_;|6jB5Z=e1H^<4}5 literal 0 HcmV?d00001 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 +