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

Asked by Bernard Ikhimwin

### Bernard Ikhimwin (view profile)

on 5 Nov 2018
Latest activity Commented on by Bernard Ikhimwin

### Bernard Ikhimwin (view profile)

on 20 Nov 2018
Accepted Answer by Image Analyst

### Image Analyst (view profile)

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.

### Image Analyst (view profile)

Answer by Image Analyst

### Image Analyst (view profile)

on 5 Nov 2018

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.

Bernard Ikhimwin

### Bernard Ikhimwin (view profile)

on 20 Nov 2018
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'
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
Image Analyst

### Image Analyst (view profile)

on 20 Nov 2018
"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
Bernard Ikhimwin

### Bernard Ikhimwin (view profile)

on 20 Nov 2018
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.