Asked by Bernard Ikhimwin
on 5 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

- I intend to look for the centre-line of an image using bwmorph operation (this I can do).
- 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).
- 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.

Answer by Image Analyst
on 5 Nov 2018

Accepted Answer

Bernard Ikhimwin
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'

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

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

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.