Detect inter-particle contact

4 views (last 30 days)
LeChat
LeChat on 1 May 2019
Hi, I have a picture of emulsion droplets. I want to detect the inter-droplet contacts and for each of these patches, plot the intensity profile.
In the following code, I first manually choose an ROI on which I try to detect one entire patch. Note that for some patches, the intensity is dimer on roughly half on the contact line: I want to measure the two plateaux intensities.
clear all
close all
fname='MTE_20190417_slice7.tif' % mind to reove the '.png' in the name of the attached file
info=imfinfo(fname);
num_images=numel(info);
flagfig=1;
k=1;% for future purposes, on a whole stack
A=imread(fname, k);
B=255-A;
figure(flagfig);imshow(B)
hold on
%% define an ROI
cursors = round(ginput(2));
plot(cursors(:,1),[1,1].*cursors(1,2),'r--')
plot(cursors(:,1),[1,1].*cursors(2,2),'r--')
plot(cursors(1,1).*[1,1],cursors(:,2),'r--')
plot(cursors(2,1).*[1,1],cursors(:,2),'r--')
%% crop the ROI
x_ROI = cursors(:,1);
y_ROI = cursors(:,2);
ROI = imcrop(B,[min(x_ROI), min(y_ROI), abs(x_ROI(2)-x_ROI(1)), abs(y_ROI(2)-y_ROI(1))]);% specifies the size and position of the crop rectangle as [xmin ymin width height]
if flagfig~=0
figure(flagfig+1)
imshow(ROI)
hold on
end
%% find the median line of the patch
binaryImage=edge(ROI,'canny');
AA=bwmorph(binaryImage, 'thicken',1);
[H,T,R] = hough(AA);
if flagfig~=0
figure(flagfig+2)
imshow(H,[],'XData',T,'YData',R,...
'InitialMagnification','fit');
xlabel('\theta'), ylabel('\rho');
axis on, axis normal, hold on;
end
P = houghpeaks(H,5,'threshold',ceil(0.8*max(H(:))));
x = T(P(:,2)); y = R(P(:,1));
if flagfig~=0
figure(flagfig+2)
plot(x,y,'s','color','white');
end
lines = houghlines(AA,T,R,P,'FillGap',5,'MinLength',7);
if flagfig~=0
figure(flagfig+3)
imshow(binaryImage), hold on
figure(flagfig+4)
imshow(AA), hold on
end
xy=nan(2,2,length(lines));
len=nan(length(lines),1);
for k = 1:length(lines)
xy(:,:,k) = [lines(k).point1; lines(k).point2];
if flagfig~=0
figure(flagfig+3)
hold on
plot(xy(:,1,k),xy(:,2,k),'LineWidth',2,'Color','green');
% Plot beginnings and ends of lines
plot(xy(1,1,k),xy(1,2,k),'x','LineWidth',2,'Color','yellow');
plot(xy(2,1,k),xy(2,2,k),'x','LineWidth',2,'Color','red');
end
% Determine the endpoints of the longest line segment
len(k) = norm(lines(k).point1 - lines(k).point2);
end
[xy_long, kmax]=max(len);% xy_long=len(kmax); % Longer length
xy_long2 =max(len(len<xy_long));% second longer length
kmax2=find(round(len*100)==round(xy_long2*100));%compute corresponding index
%
line1=xy(:,:,kmax);
line2=xy(:,:,kmax2);
linemed=(line1+line2)/2;
if flagfig~=0
figure(flagfig+1)
plot(linemed(:,1),linemed(:,2),'linewidth',2,'Color','m')
figure(flagfig+3)
plot(xy(:,1,kmax),xy(:,2,kmax),'LineWidth',2,'Color','r');
plot(xy(:,1,kmax2),xy(:,2,kmax2),'LineWidth',2,'Color','r');
end
%% compute intensity profile along this line
intprofile=improfile(255-ROI,linemed(:,1),linemed(:,2),'bilinear');
if flagfig~=0
figure(flagfig+5)
plot(intprofile)
hold on
end
I have tried several settings for the edge function, the houghpeak threshold value, the thickening of the edged image but I still have the same issue: sometimes, according to how I choose the ROI on the main image, the contact patch is not entirely detected (so the intensity profile is not complete.
Also, I still don't know how to fit the trapezoidal profile the intensity profile should (I guess) have, in order to have both the intensity values (high and low, when the latter exists).
Thank you for your help.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!