Asked by Duncan Osborne
on 2 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);

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

Geoff Hayes
on 4 Jul 2019

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]

Duncan Osborne
on 22 Jul 2019 at 11:50

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?

Geoff Hayes
on 25 Jul 2019 at 2:30

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.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## Adam (view profile)

## Direct link to this comment

https://ch.mathworks.com/matlabcentral/answers/469923-absolute-beginner-trying-to-use-a-matlab-for-a-research-internship-in-drop-shape-analysis#comment_720809

## Duncan Osborne (view profile)

## Direct link to this comment

https://ch.mathworks.com/matlabcentral/answers/469923-absolute-beginner-trying-to-use-a-matlab-for-a-research-internship-in-drop-shape-analysis#comment_720965

Sign in to comment.