From 51c9a70ad7b8b0bdfc39913eaabc6e96c6f0eb41 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Luk=C3=A1=C5=A1=20Ji=C5=99i=C5=A1t=C4=9B?= Date: Wed, 19 Jun 2024 11:26:28 +0200 Subject: [PATCH 1/1] Add the program MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit This program was created by Tereza Trávníčková. As no version control was used prior to this commit, the whole project has to be commited at once. --- .gitignore | 5 + Bubliny2023.m | 447 ++++++++++++++++++++++++++++++++++++++++++++++++++ README.txt | 20 +++ 3 files changed, 472 insertions(+) create mode 100644 .gitignore create mode 100755 Bubliny2023.m create mode 100755 README.txt diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..c1fca28 --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +*.png +*.fig +*.xls* +*.csv +Thumbs.db diff --git a/Bubliny2023.m b/Bubliny2023.m new file mode 100755 index 0000000..3c777dd --- /dev/null +++ b/Bubliny2023.m @@ -0,0 +1,447 @@ + +close all; % deletes all figures whose handles are not hidden +clear all; % removes all variables from memory +clc; % clears the Command Window display + +%% vstupni data T:\00-Orvalho\Travnickova\Bubble column +% + +%% +% User dialog for the data directory used for files + +DirectoryPath='T:\00-Orvalho\Travnickova\Bubble column\'; +DirectoryName = uigetdir(DirectoryPath,'Select the data folder'); % returns the name and path of the file selected in the dialog box +disp(['The data directory used for files is ', DirectoryName]); % displays the name of selected file in Command Window + +DDirectory = dir(DirectoryName); % Loads directories or subdirectories to struct FFiles(1).name, FFiles(1).folder, FFiles(1).date - i s casem +list={DDirectory.name}; +[Selected,OK] = listdlg('ListString',list,'ListSize',[200 500]); +SelectedDirectory = list(Selected)'; +if isempty(SelectedDirectory) % handles the error, when no spreadsheet is selected + error('No experiment selected.') +end +nSelectedDirectory = numel(SelectedDirectory); % Counts directories + +[one,Path]=size(DirectoryPath); +ResFile = [DirectoryPath DirectoryName(Path+1:end) '.xlsm']; + + +%% +% load of all directories in cycles +for iii=1:1:nSelectedDirectory % Cycles for directories with files + + +%% +% read names of all files + + FileType = '*.png'; % extension of image file + FileName=([DirectoryName '\' cell2mat(SelectedDirectory(iii)) '\' FileType]); + FFiles = dir(FileName); % Loads files to struct FFiles(1).name, FFiles(1).folder, FFiles(1).date - i s casem + nFFiles = numel(FFiles); % Counts files + +%% variable allocation + + LevelNum = zeros(1,nFFiles); + Heihgt = zeros(1,nFFiles); + Metr_fit_H = zeros(1,nFFiles); + Botom = zeros(1,nFFiles); + Metr_fit_B = zeros(1,nFFiles); + +%% +% load of all files in cycles + + step=1; + startFile= 40; %945; %200 (20. ) + endFile= 110; %nFFiles; %1000; %325 (20.) + nFFiles=endFile; % Attention!!! + for jjj=[1 startFile:step:nFFiles] % Cycles for figures files + ImageName = FFiles(jjj).name + Figure = imread([FFiles(jjj).folder '\' ImageName]); + [p,q]=size(Figure); + ContrastFigure = imadjust(Figure); % Adjust the contrast of the image using imadjust. + +%% +% Figure bands for average compution intensity, without the bands with hole and meter + + ContrastFigure1=ContrastFigure(floor(p*1/16):floor(p/8+p*1/16-3), 1:q); + [p1,q]=size(ContrastFigure1); + ContrastFigure2=ContrastFigure(floor(p/2+p/32):floor(p*3/4-p*1/16-5), 1:q); + [p2,q]=size(ContrastFigure2); + ContrastFigure3=ContrastFigure(floor(p*3/4+p*1/16):floor(p-p*1/16), 1:q); + [p3,q]=size(ContrastFigure3); + + ContrastFigureFin(1:p1,1:q)=ContrastFigure1; + ContrastFigureFin(p1+1:p1+p2,1:q)=ContrastFigure2; + ContrastFigureFin(p1+p2+1:p1+p2+p3,1:q)=ContrastFigure3; + + IntensitaAverFin{jjj}= mean(ContrastFigureFin); + + % figure + % imshow(ContrastFigure1) + % figure + % imshow(ContrastFigure2) + % figure + % imshow(ContrastFigure3) + + +%% +% Load the dimensions from the 1st image (p&q in pxl) + if jjj == 1 + +%% +% Meter + D_num=str2double(SelectedDirectory{iii}(21)); % position of camera + + reply = input('Is this a NEW EXPERIMENT that does not yet have calculated parameters for the cm scale? (y/n): ','s'); % y, n Selects where to find parameters for the cm scale. Either in the array or guess them. + if strcmp(reply,'n') + AxisIntensArray=[22,22,22,22,25,27,7,7,7,7,7,23,3,20,11,2,4,19,6,25,3]; % this needs to be guessed + AxisIntensUpArray=[60,60,65,65,34,5,29,29,29,29,29,22,26,2,27,2,28,2,29,2,29]; % this needs to be guessed + AxisEmptyArray=[190,190,190,190,200,190,190,190,190,190,190,190,190,190,190,190,190,190,190,190,190]; % this needs to be guessed 190 + LevelMeterArray=[0.95,0.95,0.95,1,0.11,0.8,0.2,0.2,0.2,0.1,0.29,0.9,0.39,0.9,0.39,0.6,0.29,0.9,0.25,0.7,0.27]; + AxisIntens=AxisIntensArray(Selected-13); % The number of all files minus 9 in which there is no data + AxisIntensUp=AxisIntensUpArray(Selected-13); + AxisEmpty=AxisEmptyArray(Selected-13); + LevelMeter=LevelMeterArray(Selected-13); + elseif strcmp(reply,'y') + AxisIntens=25; % 1.- 4. =22; 5. - 29 6.-27 7.8.9.11.11.-7 12.=23 13.=3 % Right this needs to be guessed + AxisIntensUp=2; % 3. = 65; 1. a 2. = 60; 5. - 34 6. a 7.8.9.10.11. = 29 12.=22 13.=26 % Left this needs to be guessed + AxisEmpty=190; % 5. = 200; 6. 7.8.9.10.11.12.13 = 190 % this needs to be guessed 190 + LevelMeter=0.7; % 5. = 0.11, 4. = 1; 1., 2. a 3. = 0.95, 6. =0.9; 7.8.9.=0.2 10.=0.1 11.=0.29 12.=0.9 13.=0.39 % number on the meter so that the left dash is above the surface and the right is below the surface + end + + Meter=ContrastFigure(floor(p*3/4-p*1/16+p*1/64):floor(p*3/4+p*1/64), 1:q); + ContrastMeter = imadjust(Meter); + [r,s]=size(ContrastMeter); + ContrastMeterMid=ContrastMeter(floor(end/2),:); + + hf1=figure ('Name','Contrasted meter','NumberTitle','off','PaperSize',[8.5 6],'PaperUnits','centimeters'); + imshow(ContrastMeter) + hold on + plot(ContrastMeterMid,'-m','linewidth',2) + hold on + + for kkk=2:1:q + if kkk<(1-LevelMeter)*q && ContrastMeterMid(kkk)>AxisIntensUp % looks for intensity minima for dashes on the axis above the level, Left + ContrastMeterMid(kkk)=0; + elseif kkk>=(1-LevelMeter)*q && ContrastMeterMid(kkk)>AxisIntens % looks for intensity minima for dashes on the axis under the level, Right + ContrastMeterMid(kkk)=0; + else + ContrastMeterMid(kkk)=1; + end + if ContrastMeterMid(kkk) == ContrastMeterMid(kkk-1) % looking for the right side of the peak + ContrastMeterMid(kkk-1)=0; + end + end + for mmm=q:-1:AxisEmpty + if ContrastMeterMid(mmm)==1 % erases the peaks (written numbers) to the left of the dashes on the axis + ContrastMeterMid((mmm-AxisEmpty):mmm-1)=0; + mmm=mmm-AxisEmpty; + end + end + ooo=1; + for nnn=1:1:q + if ContrastMeterMid(nnn)==1 % Writes x places of dashes on the axis + MeterScale(ooo)=nnn; + ooo=ooo+1; + end + end + + if D_num==1 + plot(MeterScale, [floor(r/2);floor(r/2);floor(r/2);floor(r/2);floor(r/2);floor(r/2);floor(r/2);floor(r/2);floor(r/2);floor(r/2)],'oc','linewidth',2) + elseif D_num==2 + plot(MeterScale(1:end-1), [floor(r/2);floor(r/2);floor(r/2);floor(r/2);floor(r/2);floor(r/2);floor(r/2);floor(r/2)],'oc','linewidth',2) + end +%-------------------------------------------------------------------------- + print ('-dpng', '-r1200', ['Contrasted_meter.png']) +%-------------------------------------------------------------------------- + + + hf2=figure ('Name','Intensity minima for dashes on the axis','NumberTitle','off','PaperSize',[8.5 6],'PaperUnits','centimeters'); + plot(ContrastMeterMid,'-m','linewidth',2) + + xlim([1,q]); + ylim([0,1]); +%-------------------------------------------------------------------------- + print ('-dpng', '-r1200', ['Intensity_places_dashes.png']) +%-------------------------------------------------------------------------- +clear Meter +clear ContrastMeter +clear ContrastMeterMid + +% Calculation of the actual position in mm using the scale of the meter + + hf3=figure ('Name','Dashes scale','NumberTitle','off','PaperSize',[8.5 6],'PaperUnits','centimeters'); + if D_num==1 + plot(MeterScale,[1000;900;800;700;600;500;400;300;200;100],'ob','linewidth',2) + hold on + polypar_Metr = polyfit(MeterScale,[1000;900;800;700;600;500;400;300;200;100],2); + elseif D_num==2 + plot(MeterScale(1:end-1),[1800;1700;1600;1500;1400;1300;1200;1100],'ob','linewidth',2) + hold on + polypar_Metr = polyfit(MeterScale(1:end-1),[1800;1700;1600;1500;1400;1300;1200;1100],2); + end + no_layer_fit=1: 1 : q; + Metr_fit=polyval(polypar_Metr,no_layer_fit); + plot(no_layer_fit,Metr_fit,'-r','linewidth',2); + hold on + if D_num==1 + xlim([0,q]); + ylim([0,1200]); + elseif D_num==2 + xlim([0,q]); + ylim([0,2000]); + end +%-------------------------------------------------------------------------- + print ('-dpng', '-r1200', 'Dashes_scale.png') +%-------------------------------------------------------------------------- + + xlswrite(ResFile,{'Contrasted meter'},['Met ' cell2mat(SelectedDirectory(iii))],'A2'); %[cell2mat(SelectedDirectory(iii)) ': Meter scale'],'A1'); ['Meter scale ' num2str(iii)] + xlswritefig(hf1, ResFile, ['Met ' cell2mat(SelectedDirectory(iii))] , 'A4') + xlswrite(ResFile,{'Intensity: places of dashes on the axis'},['Met ' cell2mat(SelectedDirectory(iii))],'A6'); + xlswritefig(hf2, ResFile, ['Met ' cell2mat(SelectedDirectory(iii))], 'B8') + xlswrite(ResFile,{'Dashes scale'},['Met ' cell2mat(SelectedDirectory(iii))],'K6'); + xlswritefig(hf3, ResFile, ['Met ' cell2mat(SelectedDirectory(iii))], 'L8') + xlswrite(ResFile,[{'Meter Scale'} {'Dashes scale'}],['Met ' cell2mat(SelectedDirectory(iii))],'A28'); + xlswrite(ResFile,[{'AxisIntens'} {'AxisIntensUp'} {'AxisEmpty'} {'LevelMeter'}],['Met ' cell2mat(SelectedDirectory(iii))],'D28'); + if D_num==1 + xlswrite(ResFile,[MeterScale' [1000;900;800;700;600;500;400;300;200;100]],['Met ' cell2mat(SelectedDirectory(iii))],'A29'); + elseif D_num==2 + xlswrite(ResFile,[MeterScale(1:end-1)' [1800;1700;1600;1500;1400;1300;1200;1100]],['Met ' cell2mat(SelectedDirectory(iii))],'A29'); + end + xlswrite(ResFile,[{AxisIntens} {AxisIntensUp} {AxisEmpty} {LevelMeter}],['Met ' cell2mat(SelectedDirectory(iii))],'D29'); + + end +%% +% Figure with average intensity, without the bands with hole and meter + figure + imshow(ContrastFigure) + hold on + plot(IntensitaAverFin{jjj},'-y','linewidth',2) + hold on + if jjj>=startFile+step %subtracts the intensities + plot(IntensitaAverFin{jjj-step},'Color',[0.9290 0.6940 0.1250],'linewidth',2) + hold on + IntensitaAverFin_2{jjj}=IntensitaAverFin{jjj}-IntensitaAverFin{jjj-step}; + waytrend=mean(IntensitaAverFin_2{jjj},2); % waytrend will tell whether to calculate the minimum or maximum green intensity IntensitaAverFin_2{jjj} + if waytrend<=0; + plot(-IntensitaAverFin_2{jjj},'-g','linewidth',2) + else + plot(IntensitaAverFin_2{jjj},'-g','linewidth',2) + end + hold on + end + + infmt = 'hh:mm:ss.SS'; + Time(jjj)=duration([(FFiles(jjj).name(8:9)) ':' (FFiles(jjj).name(10:11)) ':' (FFiles(jjj).name(12:16))],'Format',infmt)-duration([(FFiles(1).name(8:9)) ':' (FFiles(1).name(10:11)) ':' (FFiles(1).name(12:16))],'Format',infmt); + infmd = ['yy' 'MM' 'dd']; + Date(jjj) = datetime(FFiles(jjj).name(1:6),'InputFormat',infmd); + Level(jjj,1:3)=SelectedDirectory{iii}(7:9); + LevelNum(1,jjj)=str2double(Level(jjj,:))*1000; + + % % plot(lk,abs(pk),'oc','linewidth',2) + if D_num==1 + LevelPixels=MeterScale(1+(1000-LevelNum(1,jjj))/100); % x level in pixels in the original image + [pk,lk] = findpeaks(-IntensitaAverFin{jjj}(1:(LevelPixels+30))); %10 + [Heihgt(1,jjj),pkH]=max(pk); + qH=lk(pkH); + Metr_fit_H(1,jjj)=polyval(polypar_Metr,qH); + plot(qH,-Heihgt(jjj),'or','linewidth',2) + if jjj>=startFile+step + LevelPixels=MeterScale(1+(1000-LevelNum(1,jjj))/100); % x level in pixels in the original image + if waytrend<=0; + [pk,lk] = findpeaks(-IntensitaAverFin_2{jjj}(1:(LevelPixels+30))); %10 + else + [pk,lk] = findpeaks(IntensitaAverFin_2{jjj}(1:(LevelPixels+30))); %10 + end + [Heihgt(1,jjj),pkH]=max(pk); + if max(pk)>=10; % error cleaning %30 + qH=lk(pkH); + Metr_fit_H(1,jjj)=polyval(polypar_Metr,qH); + end + plot(qH,Heihgt(jjj),'om','linewidth',2) + end + + elseif D_num==2 + if jjj=startFile+step + LevelPixels=MeterScale(1); % x level in pixels in the original image + if waytrend<=0; + [pk,lk] = findpeaks(-IntensitaAverFin_2{jjj}(LevelPixels:end)); %10 + else + [pk,lk] = findpeaks(IntensitaAverFin_2{jjj}(LevelPixels:end)); %10 + end + [Heihgt(1,jjj),pkH]=max(pk); + qH=LevelPixels+lk(pkH); + Metr_fit_H(1,jjj)=polyval(polypar_Metr,qH); + if max(pk)<30; % error cleaning %30 + qH=NaN; + Metr_fit_H(1,jjj)=NaN; + end + plot(qH,Heihgt(jjj),'om','linewidth',2) + end + end + + if D_num==1 + if jjj=(Metr_fit_H(1,jjj)-2) % Throw out the bubble/liquid interface if it is above the surface, 2cm 2D interface + Metr_fit_B(1,jjj)=NaN; + qB=NaN; + elseif Botom(1,jjj)<=(max(abs(pk))+min(abs(pk)))*0.62 % 0.65 is an estimate of how much intensity the bubbles need + Metr_fit_B(1,jjj)=NaN; % For 1-4 files + qB=NaN; + end + plot(qB,Botom(jjj),'ob','linewidth',2) + end + + if jjj>=startFile+step + LevelPixels=MeterScale(1+(1000-LevelNum(1,jjj))/100); % x level in pixels in the original image + if waytrend<=0; + [pk2,lk2] = findpeaks(-IntensitaAverFin_2{jjj}(LevelPixels:end)); %10 + else + [pk2,lk2] = findpeaks(IntensitaAverFin_2{jjj}(LevelPixels:end)); %10 + end + [Botom(1,jjj),pk2B]=max(pk2); + qB=LevelPixels+lk2(pk2B); + Metr_fit_B(1,jjj)=polyval(polypar_Metr,qB); + if max(pk2)<=45; % error cleaning + qB=NaN; + Metr_fit_B(1,jjj)=NaN; + end + plot(qB,Botom(jjj),'oc','linewidth',2) + end + + elseif D_num==2 + if jjj=(Metr_fit_H(1,jjj)-2) % Throw out the bubble/liquid interface if it is above the surface, 2cm 2D interface + Metr_fit_B(1,jjj)=NaN; + qB=NaN; + elseif Botom(1,jjj)<=(max(abs(pk))+min(abs(pk)))*0.62 % 0.65 is an estimate of how much intensity the bubbles need + Metr_fit_B(1,jjj)=NaN; % For 1-4 files + qB=NaN; + end + plot(qB,Botom(jjj),'ob','linewidth',2) + end + + if jjj>=startFile+step + LevelPixels=MeterScale(1); % x level in pixels in the original image + if waytrend<=0; + [pk2,lk2] = findpeaks(-IntensitaAverFin_2{jjj}(LevelPixels:end)); %10 + else + [pk2,lk2] = findpeaks(IntensitaAverFin_2{jjj}(LevelPixels:end)); %10 + end + [Botom(1,jjj),pk2B]=max(pk2); + qB=LevelPixels+lk2(pk2B); + Metr_fit_B(1,jjj)=polyval(polypar_Metr,qB); + if max(pk2)<=20; % error cleaning + qB=NaN; + Metr_fit_B(1,jjj)=NaN; + end + plot(qB,Botom(jjj),'oc','linewidth',2) + end + end + + xlswrite(ResFile,[{'ImageName'} {'Date'} {'Time [s]'} {'Intensita Heihgt [-]'} {'Intensita Botom [-]'} {'Heihgt [cm]'} {'Botom [cm]'} {'Level [cm]'}],cell2mat(SelectedDirectory(iii)),'A1:H1'); + xlswrite(ResFile,[{ImageName} {string(Date(1,jjj))} {seconds(Time(1,jjj))} {Heihgt(1,jjj)} {Botom(1,jjj)} {Metr_fit_H(1,jjj)} {Metr_fit_B(1,jjj)} {str2double(SelectedDirectory{iii}(7:9))}], cell2mat(SelectedDirectory(iii)),['A' num2str(jjj+1) ':H' num2str(jjj+1)]); + + + clear Figure + clear ContrastFigure + clear ContrastFigure1 + clear ContrastFigure2 + clear ContrastFigure3 + clear ContrastFigureFin + % clear IntensitaAverFin + + end % jjj, Cycles for figures files + + hf4=figure ('Name','Velocity','NumberTitle','off','PaperSize',[8.5 6],'PaperUnits','centimeters'); + hp1=plot(seconds(Time(startFile:step:end)), Metr_fit_H(startFile:step:nFFiles),'-or','linewidth',2); % [1 startFile:step:nFFiles] + hold on + hp2=plot(seconds(Time(startFile:step:end)), Metr_fit_B(startFile:step:nFFiles),'-ob','linewidth',2); + hold on + hp3=plot(seconds(Time(startFile:step:end)),LevelNum(1,startFile:step:nFFiles)','-og','linewidth',2); + + xlabel('time [s]','FontWeight','bold'); + ylabel('position [cm]','FontWeight','bold'); + grid('on'); + [legend_h,object_h,plot_h,text_strings] = legend([hp1,hp2,hp3],'Heihgt', 'Botom', 'Level','Location','southeast' ); + set(legend_h, 'Box', 'off', 'FontSize',7) + +%-------------------------------------------------------------------------- + print ('-dpng', '-r1200', 'Velocity.png') +%-------------------------------------------------------------------------- + + + xlswritefig(hf4, ResFile, cell2mat(SelectedDirectory(iii)), 'J2') + +end % iii, Cycles for directories with files + diff --git a/README.txt b/README.txt new file mode 100755 index 0000000..f35e1a8 --- /dev/null +++ b/README.txt @@ -0,0 +1,20 @@ +-main program for the evaluation of experiments in BC is: +GasHoldUpSO15092023 + +two sets of data are available for evaluation (that should be possible to evaluate): +MATLAB Processing-LUK +UPW_Lab5_22062022 + +other set of data are available for evaluation (that are not possible to evaluate with the current program): +2023-08 16 - D19cm L0.8m do0.5mm + + +name of the folders with images have the following structure: +liquid phase_liquid height(in m)_superficial gas velocity(in m/s)_type of experiment(G=gas holdup; D=disengagement/shut down; E=engagement/startup)_repetition_heigh of the bottom of the column(in mm) +G.150_0.6_0.0050_G_01_0525 means: +liquid phase: 15% solution of glycerol +liquid height L = 0.6 m +superficial gas velocity q = 0.005 m/s +type of experiment: G (gas holdup) +repetition: 1st +heigh of the bottom of the column: 525 mm -- 2.30.2