MATLAB Answers

0

Diameter averaged mean intensity along the center-line of an image

Asked by Bernard Ikhimwin on 5 Nov 2018
Latest activity Commented on by Bernard Ikhimwin on 20 Nov 2018
My goal is to look for the diameter-averaged intensity of the attached image along the total center-line.
This is what I have in mind
  1. I intend to look for the centre-line of an image using bwmorph operation (this I can do).
  2. I intend to look for the branch points of the image, and subtract it from the center-line, to get discontinuous branches of the image (this I can do , using the morphological operations provided in Matlab).
  3. Finally I intend to look for the index of the pixels that are perpendicular to each pixel of the centre-line of every branch; after getting the indices of these pixels I can then use it to mask the original image, and find the mean intensity, which will be mapped to each center-line pixel of the image. I intend to plot the mean intensity as a function of the center-line.
My difficulty is in (3), I do not know how to get the index of pixels perpendicular to the center-line. I will be happy if any one can let me know how to implement this. Thank you very much as I anticipate your response.

  0 Comments

Sign in to comment.

1 Answer

Answer by Image Analyst
on 5 Nov 2018
 Accepted Answer

First get the skeleton with bwmorph(). Then get the distance transform with bwdist(). Then multiply those two images together. Then get the mean of all the non-zero elements. It should be like 3 or 4 lines of code.

  16 Comments

Dear Image Analyst, thank you for the tips, I have now coded it, and it seems to be working. However, I need more suggestions on doing somethings more efficiently. I have tried to make the code more robust, as I intend dealing with several images, Infact the previous image uploaded was just a sample image. The real ones are more complicated in terms of the interconnectedness.
Find below my code on what I have done so far, I tried it for a very simple one branched image (which I have attached). It seems to be working vaguely. My problem with complicated networks is breaking them into segments. After subtracting the branchpoints from the centerline, I still have a continuous curve, rather than segmented curve, even if it is clear that there is a junction, Matlab still recognizes it as a continuous curve. Is there a way to enforce segmentation at branch points?. Also the intensity seems very unsmooth along the centerline, and why am i getting zero intensity for some pixel, do you know how this improfile algorithm works?. Thank you
file2 = 'mysampleimage'
F1 = imread(file2,'png');
%RGB=imread (file2,'png');
I=rgb2gray(F1);
F2=im2bw(F1,0.0005);
%find the center line
F3=bwmorph(F2,'skel',Inf);
%Overlay the centerline with the bwimage
figure(1)
imshowpair(F2,F3);
%trim out any supurious branch using region of interest
F2prime=~F2+F3;
BW = roipoly(F2prime);
F3(BW)=0;
stillmoreregions = input('is there more regions you would like to mask? yes (1) and no (0) ');
while stillmoreregions>0
F2prime=~F2+F3;
BW = roipoly(F2prime);
F3(BW)=0;
stillmoreregions1= input ('is there more regions you would like to mask? yes (1) and no (0)? ');
if stillmoreregions1==0
break
end
end
%plot figure for binary image
figure(1)
imshow(F2);
%find branch points
F4=bwmorph(F3,'branchpoints');
%get line segments from centreline by subtracting branch poits from it
F5=F3-F4;
%show segments
figure(2)
imshow(F5);
%label the segments or branches
[labeledImage, numberOfObject] = bwlabel(F5,8);
figure(3)
imshow(labeledImage)
branchMeasurements = regionprops(labeledImage, 'all');
numberOfBranches = size(branchMeasurements, 1);
textFontSize = 14; % Font size
labelShiftX = -7;
for k = 1 : numberOfBranches % Loop through all branches or segments
thisBranchPixels{k} = branchMeasurements(k).PixelIdxList;
branchCentroid = branchMeasurements(k).Centroid;
fprintf(1,'#%2d', k);.
text(branchCentroid(1) + labelShiftX, branchCentroid(2), num2str(k), 'FontSize', textFontSize, 'FontWeight', 'Bold','Color','red');
end
%find the distance transform
Edist=bwdist(~F2);
figure(4)
imshow(Edist)
Rdist=Edist(F3); %radius of each pixels along the centerline
% to the boundaries of the image
%find z,y and x coordinate
x=[]; y=[]; zax=[];
for k=1:numberOfBranches
[row,col]=find(labeledImage==k);
yfic=row; xfic=col; %declare fictitious points just used to add points to x and y
x=[x xfic']; y=[y yfic'];
%get the z axis
zax1=thisBranchPixels{k};
zax=[zax zax1'];
end
figure(5)
plot(x,y)
hold on
%fit quadratic to points before or after the current point along the
% centerline, and look for a perpendicular line
windowWidth=20;
halfWindowWidth = floor(windowWidth/2);
lastIndex = length(x);
for i=1:length(x)
index1 = max(1,i - halfWindowWidth);
index2 = min(lastIndex, i + halfWindowWidth);
x1=x(index1:index2);
y1=y(index1:index2);
%get the polynomial
f = polyfit(x1,y1,2);
%find slope of line perpendicular to curve
f1=-1/(2*f(1)*x(i)+f(2));
%find right endpoints of line segment
Rxendpoint(i)=x(i)+Rdist(i)/sqrt(1+f1^2);
Ryendpoint(i)=y(i)+f1*(Rxendpoint(i)-x(i));
%find left endpoints of line sigment
Lxendpoint(i)=x(i)-Rdist(i)/sqrt(1+f1^2);
Lyendpoint(i)=y(i)+f1*(Lxendpoint(i)-x(i));
plot([Lxendpoint(i) Rxendpoint(i)],[Lyendpoint(i) Ryendpoint(i)])
end
hold off
%
%intensity along diameter
mI=zeros(length(x),1); %preallocate memory to mean Intensity mI
figure(6)
imshow(I)
Ix=double(round([Lxendpoint Rxendpoint]));
Iy=double(round([Lyendpoint Ryendpoint]));
for i=1:length(x)
q=improfile(I,Ix(2*i-1:2*i),Iy(2*i-1:2*i));
mI(i)=mean(q);
end
figure(7)
plot(zax,mI)
xlabel('zaxis')
ylabel('Intensity')
figure(10)
image(I);
colormap gray;
axis tight
axis equal
hold on
plot(Ix,Iy,'*','LineWidth',1.5)
hold off
"do you know how this improfile algorithm works?" You can edit it to see what it does. I believe it's bilinear interpolation.
edit improfile.m
Thanks Image Analyst, I was able to spot a bug in my code. I have now fixed it and it is now giving me non zero intensity. It seems to be working well. I think the only problem I will have now is processing the images well enough, so that Matlab can work with it. Thanks again.

Sign in to comment.