MATLAB Answers

0

Absolute beginner trying to use a matlab for a research internship in drop shape analysis

Asked by Duncan Osborne on 2 Jul 2019
Latest activity Commented on by Geoff Hayes
on 25 Jul 2019
Hello everyone,
I'm hoping someone very kind will help me try and figure out a piece of code I'm using for drop shape analysis. I am very sorry if my question is vague or I have missing information but I have very limited use with matlab. At the moment I am trying to resolve the problem with pure trial and error as I understand very little of the code. Thank you very much for any help.
I keep getting the following error code:
Unable to use a value of type 'ss' as an index.
Error in Example2 (line 101)
TFM.T=Rotmatrix(ACal(tilt(ss))); % generate transfer matrices using the loaded
calibration
The code is as follows:
clear all
clc
close all
% Example 2
% This example shows the processing of a video file from a tilting
% experiment. As other inputs the example needs a *.txt file with the tilt
% values of the sample (in degrees) for each frame in the video file.
% Furthermore the script loads calibration functions that describe the
% rotation and shift needed for compensating from mechanical shift of the
% stage. The calibration is loaded as three function ACal, XCal and YCal
% that are functions of tilt. Theses three functions are found by recording
% a callibration grid during tilting of the stage, points in the calbration
% grid are tracked and the shift and rotation is then calculated using the
% procrustes function without scaling. The obrained shift and rotation data
% is then fitted by cubic splines to get continious functions.
% The script can be divided into several steps
% 1. Loading all data (video, tilt values, callibration functions)
% 2. Edge detection and baseline determination in each frame
% 3. Shift baseline coordinates from each frame into coordinate system of
% first frame. Then all baseline coordinates are fitted into a single
% baseline used for all frames
% 4. Calculate contact angles for all frames using the determined baseline.
% In this part of the script there is the possibility to plot each
% frame with the corresponding contact angles and fit.
% 5. Plot the fitted contact angles and tripple line displacements as a
% function of tilt.
PlotPolynomialFrames=0;
PlotEllipticFrames=0;
%%
%---------------------------------------------------------------------------
% 1. Loading of data
%---------------------------------------------------------------------------
vidObj=VideoReader('2ul_with_annotations_jpg_without_annotations_with_compression_tif_0.avi'); %load video of tilt experiment
% load('Ex2cal.mat') %load callibration functions
fileID = fopen('2ul_with_annotations_jpg_without_annotations_with_compression_tif_0.txt');
% If you do not need the calibration you can define the calibration
% functions as:
ACal= @(t) 0*t;
XCal= @(t) 0*t;
YCal= @(t) 0*t;
% read each frame in video into variable.
mov={};
k=1;
while hasFrame(vidObj)
frame= readFrame(vidObj);
mov{k}=frame(:,:,1); %#ok<SAGROW>
k = k+1;
end
FramesTot=k-1;
% load tilt values for each frame into variable.
C = textscan(fileID,'%f');
tilt=C{1};
%%
%---------------------------------------------------------------------------
% 2. Edge detection
%---------------------------------------------------------------------------
% predefine matrices prior to loop
EdgeCell=cell(1,FramesTot); % cell structure to store detected edges
BaseVec=zeros(FramesTot,4); % matrix to store Baseline coordinates [x0L,y0L,x0R,y0R]
% run edge detection algorithm for each frame and save into variables,
% same as in exmaple 1
for n_frame=1:FramesTot
im=mov{n_frame}; %extrancting image from video
[edges, RI] = subpixelEdges(im, 25); %find edges
longestedge=findlongestedge(edges,size(im),10); % select longest edge
[edgeL,edgeR]=leftrightedges(longestedge); % divide into left and right at the drop apex
EdgeCell{n_frame,1}=edgeL; %save for later use
EdgeCell{n_frame,2}=edgeR;
[x0L,y0L,indexL]=findreflection([edgeL.x,edgeL.y],10,260); %find reflection
[x0R,y0R,indexR]=findreflection([edgeR.x,edgeR.y],10,260);
BaseVec(n_frame,:)=[x0L,y0L,x0R,y0R]; % save for later
end
%%
%---------------------------------------------------------------------------
% 3. Average baselines found in all frames.
%---------------------------------------------------------------------------
% Define functions used to perform coordinate rotation
% Z being the coordinate system of the first frame
% Y being local coordinate system in each frame.
Rotmatrix=@(angle) [cosd(angle),-sind(angle);sind(angle),cosd(angle)];
Coordinates2Z=@(coordinates,TransferMatrix) bsxfun(@plus,coordinates*TransferMatrix.T,TransferMatrix.c);
Coordinates2Y=@(coordinates,TransferMatrix) bsxfun(@minus,coordinates,TransferMatrix.c)/(TransferMatrix.T);
for BaseVecZ=zeros(size(BaseVec));
TFM.T=Rotmatrix(ACal(tilt(ss))); % generate transfer matrices using the loaded calibration
TFM.c=[XCal(tilt(ss)),YCal(tilt(ss))];
BaseVecZ(ss,1:2)=Coordinates2Z(BaseVec(ss,1:2),TFM); % Shift coordinate system of baselines to coordinate system of first frame (Z)
BaseVecZ(ss,3:4)=Coordinates2Z(BaseVec(ss,3:4),TFM);
end
% Find the average baseline using all baseline coordinates
xBase=[BaseVecZ(:,1);BaseVecZ(:,3)];
yBase=[BaseVecZ(:,2);BaseVecZ(:,4)];
B=AverageBaseline(xBase,yBase);
AverageBaseVecZ=bsxfun(@times,ones(FramesTot,4),[0,B(1),size(im,2),size(im,2)*B(2)+B(1)]);
% The new baseline is found in the Z coordinate system and is transfered to
% the Y coordinate system of each frame to be used for contact angle
% measurements.
AverageBaseVecY=zeros(FramesTot,4);
for ss=1:FramesTot
TFM.T=Rotmatrix(ACal(tilt(ss)));
TFM.c=[XCal(tilt(ss)),YCal(tilt(ss))];
AverageBaseVecY(ss,1:2)=Coordinates2Y(AverageBaseVecZ(ss,1:2),TFM);
AverageBaseVecY(ss,3:4)=Coordinates2Y(AverageBaseVecZ(ss,3:4),TFM);
end
%%
%---------------------------------------------------------------------------
% 4. Fit contact angles using averaged baseline
%---------------------------------------------------------------------------
% predefine variables prior to loop
TLvec=zeros(FramesTot,4);
PolyCAS=zeros(FramesTot,2);
EllipseCAS=zeros(FramesTot,2);
% fit contact angles using the edges stored in EdgeCell as in example 1
for n_frame=1:FramesTot
edgeL=EdgeCell{n_frame,1};
edgeR=EdgeCell{n_frame,2};
v=num2cell(AverageBaseVecY(n_frame,:));
[x0L,y0L,x0R,y0R]=v{:};
PolyData=polynomialfit(edgeL,edgeR,[x0L,y0L],[x0R,y0R]);
EllipseData=EllipticFit(edgeL,edgeR,[x0L,y0L],[x0R,y0R]);
% Store obtained data for later use
PolyCAS(n_frame,1)=PolyData.CAL;
PolyCAS(n_frame,2)=PolyData.CAR;
EllipseCAS(n_frame,1)=EllipseData.CAL;
EllipseCAS(n_frame,2)=EllipseData.CAR;
TLvec(n_frame,1:2)=EllipseData.TLL;
TLvec(n_frame,3:4)=EllipseData.TLR;
% Option to plot
if PlotPolynomialFrames % plot polynomial data
figure %#ok<UNRCH>
imshow(mov{n_frame})
hold on
t=linspace(-3,3);
plot((x0R-x0L)*t+x0L,(y0R-y0L)*t+y0L,'r--','LineWidth',2)
radius=60;
imrotation=atand((y0R-y0L)/(x0R-x0L));
lw=2;
plot([PolyData.TLR(1),PolyData.TLR(1)-radius*cosd(PolyData.CAR+imrotation)],[PolyData.TLR(2),PolyData.TLR(2)-radius*sind(PolyData.CAR+imrotation)],'g--','LineWidth', lw)
plot([PolyData.TLL(1),PolyData.TLL(1)+radius*cosd(PolyData.CAL-imrotation)],[PolyData.TLL(2),PolyData.TLL(2)-radius*sind(PolyData.CAL-imrotation)],'g--','LineWidth', lw)
plot(PolyData.EvalPolyR(:,1),PolyData.EvalPolyR(:,2),'g--','LineWidth',2)
plot(PolyData.EvalPolyL(:,1),PolyData.EvalPolyL(:,2),'g--','LineWidth',2)
end
if PlotEllipticFrames % plot ellipse data
figure; %#ok<UNRCH>
imshow(mov{n_frame})
hold on
t=linspace(-3,3);
plot((x0R-x0L)*t+x0L,(y0R-y0L)*t+y0L,'r--','LineWidth',2)
radius=60;
imrotation=atand((y0R-y0L)/(x0R-x0L));
ellipserim=@(e,t) ones(size(t))*[e.X0_in,e.Y0_in]+cos(t)*[cos(-e.phi),sin(-e.phi)]*e.a+sin(t)*[-sin(-e.phi),cos(-e.phi)]*e.b;
lw=2;
plot([EllipseData.TLR(1),EllipseData.TLR(1)-radius*cosd(EllipseData.CAR+imrotation)],[EllipseData.TLR(2),EllipseData.TLR(2)-radius*sind(EllipseData.CAR+imrotation)],'b-','LineWidth', lw)
plot([EllipseData.TLL(1),EllipseData.TLL(1)+radius*cosd(EllipseData.CAL-imrotation)],[EllipseData.TLL(2),EllipseData.TLL(2)-radius*sind(EllipseData.CAL-imrotation)],'b-','LineWidth', lw)
PlotEllipseDataL=ellipserim(EllipseData.ellipse{1},linspace(0,2*pi,1000)');
PlotEllipseDataR=ellipserim(EllipseData.ellipse{2},linspace(0,2*pi,1000)');
plot(PlotEllipseDataL(:,1),PlotEllipseDataL(:,2),'b','LineWidth',2)
plot(PlotEllipseDataR(:,1),PlotEllipseDataR(:,2),'b','LineWidth',2)
end
end
% Caluclate the tripple line position in Z coordinate system and calculate
% displacement from first frames as a function of tilt.
shifted_base_vec=zeros(size(TLvec));
for ss=1:FramesTot
TFM.T=Rotmatrix(ACal(tilt(ss)));
TFM.c=[XCal(tilt(ss)),YCal(tilt(ss))];
shifted_base_vec(ss,1:2)=Coordinates2Z(TLvec(ss,1:2),TFM);
shifted_base_vec(ss,3:4)=Coordinates2Z(TLvec(ss,3:4),TFM);
end
baseshift=bsxfun(@minus,shifted_base_vec,shifted_base_vec(1,:));
distance_vec=[sqrt(baseshift(:,1).^2+baseshift(:,2).^2),sqrt(baseshift(:,3).^2+baseshift(:,4).^2)];
%%
%---------------------------------------------------------------------------
% 5. Plot result for video analysis
%---------------------------------------------------------------------------
% the figure is divided into two a lower (ax1) figure with the contac angles and
% an upper (ax2) with the displacement of the tripple line.
scale=2;
FigHandle=figure;
set(FigHandle,'Units','Centimeters');
set(FigHandle,'Position',[4, 1, 8*scale, 6*scale]);
ax1=axes('Units','Centimeters');
hold on
p1=plot(tilt,PolyCAS(:,1),'rx--');
p2=plot(tilt,PolyCAS(:,2),'rx-');
p3=plot(tilt,EllipseCAS(:,1),'bo--');
p4=plot(tilt,EllipseCAS(:,2),'bo-');
hYLabel1=ylabel('Contact angle [degrees]');
hXLabel=xlabel('Tilt [degrees]');
axespos=get(ax1,'Position');
set(FigHandle,'Position',[4, 1, 8*scale, 6*scale*1.2]);
set(ax1,'Position',[axespos(1:3),axespos(4)*0.8])
ax2=axes;
set(ax2,'Units','Centimeters','Position',[axespos(1),axespos(2)+axespos(4)*0.8*1.05,axespos(3),axespos(4)*0.8*0.5]);
hold on
p5=plot(tilt,distance_vec(:,1),'k--');
p6=plot(tilt,distance_vec(:,2),'k-');
% olny for legend purposes
p7=plot(tilt,zeros(size(tilt))-10,'--k');
p8=plot(tilt,zeros(size(tilt))-10,'-k');
p9=plot(tilt,zeros(size(tilt))-10,'rx');
p10=plot(tilt,zeros(size(tilt))-10,'bo');
ax2.YLim=[0 ceil(max(max(distance_vec))/10)*10];
hYLabel2=ylabel('Tripple line displacement [pixel]');
set(ax2,'XaxisLocation','top')
ax2.XTickLabel=cell(size(ax2.XTickLabel,1),1,'');
set([ax1,ax2],'Box','off')
hLegend=legend([p7,p8,p9,p10],{'left','right','Polynomial','Ellipse'},'Location','NW');
ax1.XLim=[0 ceil(max(max(tilt)))];
ax2.XLim=[0 ceil(max(max(tilt)))];
% The rest is cosmetic changes only
set([ax1,ax2],'FontName','Helvetica');
t1=get(ax1,'children');
set(t1,'LineWidth',0.5*scale)
t2=get(ax2,'children');
set(t2,'LineWidth',0.5*scale)
set([ax1,ax2],...
'Box','off',...
'TickDir','out',...
'TickLength',[.02 .02],...
'XMinorTick','on',...
'YMinorTick','on',...
'YGrid','on',...
'XColor',[.3 .3 .3],...
'YColor',[.3 .3 .3],...
'LineWidth',1*scale);

  2 Comments

It's very useful, as a beginner, to familiarise yourself with debugging options. They are very simple to use in Matlab, in general, but powerful for finding these bugs and being able to check things on command line.
In particular, the 'Pause on errors' option from the breakpoints menu (or equivalent code version) is exceptionally useful to stop the code when an error occurs and allow you to investigate whilst still in the code itself at that point rather than the usual behaviour where an error just stops code execution and dumps you out of any functions that had the error.

Sign in to comment.

2 Answers

Answer by Geoff Hayes
on 2 Jul 2019
 Accepted Answer

Duncan - this block of code seems suspect
for BaseVecZ=zeros(size(BaseVec));
TFM.T=Rotmatrix(ACal(tilt(ss))); % generate transfer matrices using the loaded calibration
TFM.c=[XCal(tilt(ss)),YCal(tilt(ss))];
BaseVecZ(ss,1:2)=Coordinates2Z(BaseVec(ss,1:2),TFM); % Shift coordinate system of baselines to coordinate system of first frame (Z)
BaseVecZ(ss,3:4)=Coordinates2Z(BaseVec(ss,3:4),TFM);
end
ss hasn't been defined prior to this block of code, and the assignment in the for line does not make sense. Perhaps this should instead be
BaseVecZ=zeros(size(BaseVec));
for ss = 1:size(BaseVec, 1)
TFM.T=Rotmatrix(ACal(tilt(ss))); % generate transfer matrices using the loaded calibration
TFM.c=[XCal(tilt(ss)),YCal(tilt(ss))];
BaseVecZ(ss,1:2)=Coordinates2Z(BaseVec(ss,1:2),TFM); % Shift coordinate system of baselines to coordinate system of first frame (Z)
BaseVecZ(ss,3:4)=Coordinates2Z(BaseVec(ss,3:4),TFM);
end

  16 Comments

try creating the file from the MATLAB file editor. And try the attached file too. I created this file through MATLAB and was able to read the three values from it.
>> fileID2 = fopen('geoff.txt');
>> D = textscan(fileID2,'%f');
>> D
D =
[3x1 double]
Hi Geoff,
Apologies for the delay I've been away. Using the other text file resolved the specific issue, however I'm not getting another.
Error using textscan
Invalid file identifier. Use fopen to generate a valid file identifier.
Error in Example2 (line 62)
C = textscan(fileID,'%f');
I'm assuming it doesn't like the type of file I'm using. I've had a play around with different types of file but with no luck. Is there anything you could reccomend?
Duncan - given the error message, the problem could be that the file cannot be found. Try including the full path (to the file) along with the file name. i.e.
fileID2 = fopen('/users/geoffh/documents/geoff.txt');

Sign in to comment.


Answer by Image Analyst
on 4 Jul 2019

tilt does not seem to be a vector so you can't use ss to index into it.

  0 Comments

Sign in to comment.