%%% Determine interface thickness from shadowgraph still images %%% %%% This method is manual and is not extremely accurate and entirely %%% repeatible but it gives valuable information. %%% %%% Adrien Lefauve, February 2018 clear all close all % Import jpg image [filename, user_canceled] = imgetfile('MultiSelect',1); %select the image(s) to analyse n_im = size(filename,2); for i_im=1:n_im [img,cmap] = imread(filename{i_im}); if size(img,3) == 1 img=ind2gray(img,cmap); elseif size(img,3) == 3 img=rgb2gray(img); end imgs(:,:,i_im) = double(img)/255; end im = squeeze(mean(imgs,3)); % Determine angle of duct to rotate the image back to horizontal scrsz = get(0, 'ScreenSize'); set(gcf,'Position',scrsz*0.9); imshow(im); % colormap('bone'); %show image title(['Click 2 times along top duct boundary (left and right hand sides) to determine angle, and 1 time along bottom boundary (anywhere)'],'FontSize',18); [x, z] = ginput(3); close all fit = polyfit(x(1:2),z(1:2),1); phi = atan(fit(1))*180/pi; %angle in degrees if( abs(phi)>1 ) %as this will mess up image slightly, only do if necessary (angle is >1 degree) im_rot = imrotate(im,phi,'bilinear'); else im_rot = im; end [x_px,z_px] = size(im_rot); z_max = floor(z(3))+50; z_min = floor(min(z(1:2)))-50; % Chose 1 location (where profile will be visible, i.e. no big bubble in the way) figure scrsz = get(0, 'ScreenSize'); set(gcf,'Position',[1950 200 1600 800]); imshow(im_rot); %colormap(cmap); %show image hold on; title('Select 1 location by clicking','FontSize',14); [x_pos,~] = ginput(1); x_pos=floor(x_pos); close all % Determine thickness at this location hCurr = figure scrsz = get(0, 'ScreenSize'); set(gcf,'Position',scrsz*0.9); imshow(im_rot(z_min:z_max,50:end-400,:)); hold on; plot([x_pos,x_pos],[1,z_max-z_min+1],'g','LineWidth',2); avg = mean(mean(im)); intensity = im_rot(z_min:z_max,x_pos,1); %NB: we're only reading one colour intensity... but that shouldn't matter too much really! intensity = smooth(intensity,3); avg = mean(intensity); intensity = intensity - avg; intensity = intensity / max(abs(intensity)); plot(x_pos + 100*intensity,1:z_max-z_min+1,'c'); title(['Click on top boundary, top layer, bottom layer, bottom boudary'],'FontSize',14); [x, z] = ginput(4); close; z_all = z + z_min - 1; close all relative_thickness = (z_all(2,:)-z_all(3,:))./(z_all(1,:)-z_all(4,:)) % in percents of full duct height H (%) % This data is then copied to an excel spreadsheet and to a .mat structure % (also supplied in the data repository) % The code below was used to plot the delta(theta,Re) data as in the % paper, and to plot the least-squares fit as black contours % % %% ADD CONTOUR OF FIT % % mSID % p00 = 2.459; % p10 = 2.545; % p01 = -0.5315; % p20 = 0.416; % p11 = -0.2042; % p02 = 0.1393; % % % tSID % p00 = 2.122; % p10 = -0.6107; % p01 = -1.58; % p20 = 0.5708; % p11 = 0.7026; % p02 = 0.4047; % % % LSID % p00 = -2.02 ; % p10 = 2.51 ; % p01 = 2.254 ; % p20 = 0.274 ; % p11 = -0.2708 ; % p02 = -0.3109 ; % % % HSID % p00 = 1.877 ; % p10 = 2.375 ; % p01 = -0.0453 ; % p20 = 0.145 ; % p11 = -0.3149 ; % p02 = 0.005809 ; % % [X,Y]=meshgrid(-2:0.01:-1,2.9:0.01:4.5); % f= p00 + p10*X+ p01*Y + p20*X.^2 + p11*X.*Y + p02*Y.^2; % % figure % opt.Markers = {}; % markersize = 30; % hold on % [C,h]=contourf(10.^X,10.^Y,f,[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9],'-k','Showtext','off','Linewidth',1); caxis([0,1]); colormap(cmap_bw); % % clabel(C,h,'FontSize',12,'Color','k','FontName','Times') % % % % SCATTER PLOT OF DATA ON TOP % cmap_bw = kron([1 1 1],linspace(1,0,100)'); % p1 = scatter(sind(theta),Re,markersize,dI,'filled'); colormap(cmap_bw); caxis([0,1]); % hold on % scatter(sind(theta),Re,markersize,'ok'); % opt.Markers = [opt.Markers ,'o']; % % % % Settings to make plot nicer % opt = []; % opt.BoxDim = [5 3.8]; % opt.FontName = 'Times'; % % opt.XLabel = '$\sin \theta $'; % xlabel % opt.XScale = 'log'; % opt.XLim = [1e-3 2.5e-1]; % opt.XTick = [1e-3 1e-2 1e-1]; % opt.XGrid = 'on'; % opt.XMinorTick = 'off'; % % opt.YLabel = '$Re$'; %ylabel % opt.YScale = 'log'; % opt.YTick = [1e2, 1e3,1e4,1e5]; % opt.YLim = [1e2 1e5]; % opt.YGrid = 'on'; % % % opt.LegendBox = 'off'; % % if iind == 3 % % opt.Legend = {'H', 'I', 'T'}; % % elseif iind == 4 % % opt.Legend = {'L', 'H', 'I', 'T'}; % % elseif iind == 5 % % opt.Legend = {'L', 'W', 'H', 'I', 'T'}; % % end % % opt.LegendBoxColor = 'w'; % % opt.LegendLoc = 'NorthEast'; % opt.AxisLineWidth = 1; % % % % % Make plot nicer % setPlotProp(opt); % % % Adjust Y label % set(get(gca,'YLabel'),'Rotation',0) % set(gca,'Position',get(gca,'Position')+[0.3 0.1 0 0]) % % xh = get(gca,'Ylabel'); % handle to the label object % p = get(xh,'position'); % get the current position property % p(1) = 0.75*p(1) ; % p(2) = 1.1*p(2) ; % % negative values put the label below the axis % set(xh,'position',p) % set the new position % % xh = get(gca,'Xlabel'); % handle to the label object % p = get(xh,'position'); % get the current position property % p(2) = 0.9*p(2) ; % double the distance, % % negative values put the label below the axis % set(xh,'position',p) % set the new position % % % % Save it % set(gcf,'Units','inches'); % screenposition = get(gcf,'Position'); % set(gcf,... % 'PaperPosition',[0.1 0.1 screenposition(3:4)],... % 'PaperSize',[screenposition(3:4)+0.1]); % % % print -dpdf flux % movefile('flux.pdf', strcat(char('interface-thickness'),'.pdf'));