Calculate the Cone angle of a spray image

Hi, I'm having some problems in calculating the cone angle for an array of diesel sprays. i tried initially using the bwtraceboundary function, however, this has problems in analysing multiple images. At the very start of the spray development it does not work(when the spray is very small) and in the late stages (when the injector has finished spraying, the spray is large, in contact with the far wall and no longer in contact with the nozzle)it also does not work.

Below is one of the original images and its cropped and thresholded version.

after re-sizing and thresholding the image, the matlab code i used was as follows:

dim=size(BW);

%left edge

col1=56;

row1=min(find(BW(:,col1), 1 ));

%right edge

col2=63;

row2=min(find(BW(col2,:), 1 ));

boundary1=bwtraceboundary(BW,[row1,col1],'N',8,220,'counter');

boundary2=bwtraceboundary(BW,[row2,col2],'N',8,220);

% fit points to line

ab1=polyfit(boundary1(:,2),boundary1(:,1),1);

ab2=polyfit(boundary2(:,2),boundary2(:,1),1);

r = ab1(1) .* boundary1(:,2) + ab1(2)+offsetY;

% now plot both the points in y and the curve fit in r onto original image using an offset

hold on;

plot(offsetX+boundary1(:,2), r, 'g','LineWidth',2);

r = ab2(1) .* boundary2(:,2) + ab2(2)+offsetY;

hold on;

plot(offsetX+boundary2(:,2), r, 'g','LineWidth',2);

% create vectors

vect1=[1 ab1(1)];

vect2=[1 ab2(1)];

dp=dot(vect1,vect2);

length1=sqrt(sum(vect1.^2));

length2=sqrt(sum(vect2.^2));

% obtain angle

angle=180-acos(dp/(length1*length2))*180/pi;

%This was then plotted on top of the original sized image

However, I am sure there are many other far more accurate and consistent ways of doing this. My knowledge of matlab is very poor, i have used tutorials up until now, but there don't seem to be any more that can help. Although I know that a good way of doing this is by scanning the image in a certain direction and finding the boundary between white and black pixels, unfortunately i do not know how to write any of this in code.

Any help would be really appreciated!

Thank you very much

Jed

[지영 허]

Hello Jed

I'm so happy to find the code you made. I almost succeeded in getting the spray cone angle through the code you made. However, it is difficult to complete, so I would like to practice with your picture and code, so could you tell me all of your codes? I hope there is a part like offsetX or offsetY missing. I'd really appreciate it if you could let me know.

Answer 1

[Image Analyst]

No, you don't want to use bwtraceboundary. You can use bwboundaries instead if you can isolate the plume and get rid of the surrounding circle so you have a nice binary mask of just the plume/spray.

Assuming you've identified the apex of the triangle at the top, then the spray angle changes depending on how far down you want to measure it. And it will be different if you just look at a line and use that width, or if you do a fit from all the side points from the current line back up to the origin point. So do you want the "instantaneous" width plotted at a function of line down from the origin, or do you want the "fitted width"? If you go for the instantaneous width, then I guess you could just get the widest angle. If you get the fit, then you need to decide where to stop the fit because it starts to taper down again near the bottom.

[Jed]

I have got an apex which will be specified by the user by entering values for InjectorX & InjectorY to give a coordinate to be applied to all of their images. I ideally wanted a line of best fit of the points between the origin and the point of maximum width on either side which i have tried to show in the attempted version above, but i understand that is very inconsistent as you have to change the line length accordingly.I just really wanted something that would show consistency between multiple images for good comparison but with my very basic knowledge this has been very tricky.

[Image Analyst]

Makes sense and doesn't seem very hard. Just scan down line by line and find the left most and right most pixel on each line. Subtract them and get the width. Note the line that has the widest width. Then just send the x and y values into polyfit(). Do for each side to get the fitted line for each side, and then do some trig to get the angle between the lines. You should be able to do this.

[Jed]

Thank you very much! I think i will be alright doing the polyfit and the angle calculation but it is the coding of the scan which i do not know how to do. Would i be able to use your previous code of:

for columnNumber = 1:size(ImageArray, 2)

    topRows(columnNumber) = find(ImageArray(:, columnNumber)>0, 1,'first');

    bottomRows(columnNumber) = find(ImageArray(:, columnNumber)>0, 1,'last');

end

and just use the part for topRows due to my orientation? or would this not work due to the surrounding circle?

[Image Analyst]

Yes but you'd have to rotate your image 90 degrees. If you want to just crop it and leave the plume vertical then adjust your code to go down rows:

for row = 1:size(binaryImage, 1)
    leftColumns(row) = find(binaryImage(row, :)>0, 1,'first');
    rightColumns(row) = find(binaryImage(row, :)>0, 1,'last');
end

[Jed]

Aw ok, thanks. I have tried it both rotated and vertical using the following code(altering it for the orientation and adding an offset in both x and y so i can impose it on top of the original image):

ab1=polyfit(leftColumns(:,2),leftColumns(:,1),1);
ab2=polyfit(rightColumns(:,2),rightColumns(:,1),1);

%plot leftside
r = ab1(1) .* leftColumns(:,2) + ab1(2)+offsetY; 
hold on;
plot(offsetX+leftColumns(:,2), r, 'g','LineWidth',2);

%plot rightside
r = ab2(1) .* rightColumns(:,2) + ab2(2)+offsetY;
hold on;
plot(offsetX+rightColumns(:,2), r, 'g','LineWidth',2);

but when trying either i get the error message:

Warning: Polynomial is not unique; degree >= number of data points. 
> In polyfit at 71
  In MeasureAngleForum at 59 
Warning: Polynomial is not unique; degree >= number of data points. 
> In polyfit at 71
  In MeasureAngleForum at 60

Is this due to the way i have tried to implement polyfit?

[Image Analyst]

Maybe try to have your "x" direction be the direction along the plume and the "y" direction be the cross direction. Maybe it doesn't like having 2 Y for the same x. If you want, attach a script that gets the binary image of the plume and then does your code for finding edges and fitting.

[Jed]

Thank you very much for all of your time. Yea i get the feeling it is because it doesn't like having two Y values. I have no idea why though.

My entire script is as follows:

% load image
cd('H:\Individual Project\MATLAB\SprayExamples\2');
RG=imread('2011.jpg');
RGB=imrotate(RG,271.6);
imshow(RGB)

% crop image
start_row=210;
start_col=170;
cropRGB=RGB(start_row:310,start_col:440);

offsetX=start_col-1;
offsetY=start_row-1;

% threshold
t=200;
m=size(cropRGB,1); 
n=size(cropRGB,2);

BW=zeros(m,n); 
for y=1:n;
for x=1:m;
if cropRGB(x,y)>t
BW(x,y)=255; 
else
BW(x,y)=0; 
end
end
end
imshow(BW);

%find boundaries of spray
for columnNumber = 1:size(BW, 2)
  topRows(columnNumber) = find(BW(:, columnNumber)>0, 1,'first');
  bottomRows(columnNumber) = find(BW(:, columnNumber)>0, 1,'last');
end
 ab1=polyfit(topRows(:,2),topRows(:,1),1);
 ab2=polyfit(bottomRows(:,2),bottomRows(:,1),1);

%impose on top of thresholded, rotated original
% threshold original
t=180;
m=size(RGB,1); 
n=size(RGB,2);

BW2=zeros(m,n); 
for y=1:n;
for x=1:m;
if RGB(x,y)>t
BW2(x,y)=255; 
else
BW2(x,y)=0; 
end
end
end
imshow(BW2);
hold on;

%plot topside
r = ab1(1) .* topRows(:,2) + ab1(2)+offsetY; 
hold on;
plot(offsetX+topRows(:,2), r, 'g','LineWidth',2);

 %plot bottomside
  r = ab2(1) .* bottomRows(:,2) + ab2(2)+offsetY;
  hold on;
  plot(offsetX+bottomRows(:,2), r, 'g','LineWidth',2);

But afterwards it just displays that same error message...

Answer 2

[Sean de Wolski]

Here's the approach I would take. Segment the plume, sum up the rows, calculate the slope of the rows up to the peaks, polyfit that.

With your original image being '2011.jpg'

%%Load, green slice

I = imread('2011.jpg');
green = I(:,:,2);

%%Find plume

P = (im2bw(green,graythresh(green)));
P = xor(P,imfill(P,'holes'));
P = keepMaxObj(P);

% imshow(P)

%%Engine

npx = sum(P,2); % sum along second dimension
npx = nonzeros(npx); % extract nonzeros
[~,idx] = max(npx);
npx = npx(1:idx);
plot(npx);
mb = polyfit((1:idx).',npx,1); % fit line

slope = mb(1);
the_angle = atand(slope)

[Sean de Wolski]

Oh, and here's keepMaxObj.

[Jed]

Thank you very much!keepMaxObj is brilliant. I'm afraid i don't follow what goes on after %%Engine though, as it produces a graph and a single value for slope?

[Jed]

Also, i wouldn't be able to use keepMaxObj on later images where the spray makes contact with the surrounding circle, so would i be best off cropping all of the images first so this can't happen?

[Sean de Wolski]

One thing, you don't actually need to halve mb(1). mb(1) would be the slope relative to one side which is the slope. that should be fed into atand() see fix.

What happens is I take the number of pixels in each row that are part of the plume, i.e. the width of the plume, and then fit a line to that using polyfit. The plot was just to visualize the thickness of each row up to the maximum thickness.

[Sean de Wolski]

As for generalization, it's hard to say without seeing another image.

[Jed]

Aw ok, that makes more sense now! thank you! The images range from having a nonexistent spray to having a fully developed spray where the plume is not in contact with the nozzle. Is there a way of applying this to images like the ones attached below? When i've tried it with other images it just produces a slope of 1 and angle of 45deg with no plume detected. Also, is it possible to plot lines showing the tangent to spray/angle measured on top of the original image for presentational purposes?

Answer 3

[bimal khatiwada]

Hi everyone,

Can you please tell me the MATLAB code to calculate spray length and spray cone angle for any images? For example:

[Image Analyst]

Didn't Sean's code do that?

https://www.mathworks.com/matlabcentral/answers/120866-calculate-the-cone-angle-of-a-spray-image#answer_127990