2 views (last 30 days)

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);

Geoff Hayes
on 2 Jul 2019

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]

Geoff Hayes
on 25 Jul 2019

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.

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.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## Direct link to this comment

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

⋮## Direct link to this comment

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

## Direct link to this comment

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

⋮## Direct link to this comment

https://se.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.