Calculate the radial pixeldistance

Hi all,

I have an image which contains spots. I want to calculate radial pixeldistance from the surface and 20 pixels from it. Fact is that the grayscale value will change as we move away from the "spot".

or

calculate the circumferential mineral density (from the "surface" of an osteocyte) as a function of the radial pixeldistance from the osteocytes and up to 20 pixels from the osteocytes.

Does anyone know a solution to this "problem".

/Best regards

Answer 1

[Image Analyst]

I've done something similar where I had to find out the dye penetration into the spots as a function of distance into the spots, and give an average dye penetration profile - a 2D plot of basically gray level as a function of distance from the outer edge.

You need to compute the Euclidean Distance Transform with bwdist(). It sounds like you want the gray levels outside the white spots out to a distance of 20 pixels. So take the EDT and threshold it at 20. Now you have an original gray scale image, and an EDT image that says how far the pixel is from any white spot. Now you can just scan the pixels in the EDT image and create a histogram of the gray levels for every distance. Pretty easy. Give it a try.

[mehdi]

we are doing this below code.

D = bwdist(BW);

figure; str = datestr(now,0);

[c, h]=contourf (D,8);

colormap jet; alpha (.3);

title([' Generated on: ',str]);

xlabel('x-position');

ylabel ('y-position');

clabel(c,h);

how you mean we can find distance of 20 pixel for each spots.

we are finding value for each pixels around spots. but how to find distance to 20th pixels??

[Image Analyst]

I don't understand. Are you working together with Kenni? Or doing the same thing he wants to do?

I can't really do much since Kenni didn't upload the original gray scale image, just the segmented, binary image.

[mehdi]

yes we are working toghether. we are in same group that making a project.

here is our whole code:

[Filename, Pathname] = uigetfile(... {'*.bmp;*.jpg;*.tif;*.png;*.gif','Image files (.bmp;.jpg;*.tif;*.png;*.gif)';... '*.*','All file (.)'},... 'Select a file to open');

if ~isequal(Filename,0)

FileName = fullfile(Pathname,Filename);

else

 error('User pressed cancel.');

end

im = imread(FileName);

figure (1) imagesc(im);

I = imcrop; imagesc(I);

level = graythresh(I);

BW = ~im2bw (I, level);

BW = imclearborder(BW);

BW = bwareaopen(BW,50);

BW = imfill(BW,'holes');

L = bwlabeln(BW);

s = regionprops(L,'Area');

Areal=[s.Area];

idx = find(Areal >= 500);

BW(ismember(L, idx)) = 0;

figure (2) imagesc(BW);

D = bwdist(BW);

figure(3);

str = datestr(now,0);

[c, h]=contourf (D,20);

colormap jet; alpha (.3);

title([' Generated on: ',str]);

xlabel('x-position');

ylabel ('y-position');

clabel(c,h);

we hope you can help us with it.

thanks.

[Image Analyst]

Well I was going to but I see you didn't attach your image like I asked. So I'll check back later tonight perhaps, after I return. Or just specify one of the standard demo images that I can use, like rice.png or something.

[Image Analyst]

OK, I created a sample image using the peaks function, then thresholded it to create a binary image. Then I did what I said to get the average radial profile, which is the mean gray level as a function of distance from the nearest white blob. The code below produces the image below.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

% Create some sample data.
[X,Y,grayImage] = peaks(330);
grayImage = uint8(255 * mat2gray(grayImage));
% Display the original image.
subplot(2, 3, 1);
imshow(grayImage, []);
axis on;
title('Original Image', 'FontSize', fontSize);
set(gcf, 'units','normalized','outerposition',[0 0 1 1]); % Maximize figure.

% Let's compute and display the histogram.
[pixelCount, grayLevels] = imhist(grayImage);
subplot(2, 3, 2); 
bar(grayLevels, pixelCount);
grid on;
title('Histogram of Original Image', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.

% Binarize it to find bright things, such as the dike.
binaryImage = grayImage > 160;
% Display the image.
subplot(2, 3, 3);
imshow(binaryImage, []);
title('Binarized', 'FontSize', fontSize);
drawnow;

% Compute the Euclidean Distance Transform
edtImage = bwdist(binaryImage);
% Display the image.
subplot(2, 3, 4);
imshow(edtImage, []);
title('Euclidean Distance Transform', 'FontSize', fontSize);

% Get the histogram of it so we know what distances are there
% allDistances = unique(edtImage)
maxDistance = ceil(max(edtImage(:)))

% Let's compute and display the histogram.
[pixelCounts, binCenters] = hist(edtImage(:), 100);
subplot(2, 3, 5); 
bar(binCenters, pixelCounts);
grid on;
title('Histogram of EDT image', 'FontSize', fontSize);
xlim([0 binCenters(end)]); % Scale x axis manually.

% Allocate an array for the profile
profileSums = zeros(1, maxDistance);
profileCounts = zeros(1, maxDistance);
% Scan the original image getting gray level, and scan edtImage getting distance.
% Then add those values to the profile.
[rows, columns] = size(edtImage);
for column = 1 : columns
  for row = 1 : rows
    thisDistance = round(edtImage(row, column));
    if thisDistance <= 0
      continue;
    end
    profileSums(thisDistance) = profileSums(thisDistance) + double(grayImage(row, column));
    profileCounts(thisDistance) = profileCounts(thisDistance) + 1;
  end
end
% Divide the sums by the counts at each distance to get the average profile
averageRadialProfile = profileSums ./ profileCounts;
subplot(2, 3, 6); 
plot(1:length(averageRadialProfile), averageRadialProfile, 'b-', 'LineWidth', 3);
grid on;
title('Average Radial Profile', 'FontSize', fontSize);
xlabel('Distance from nearest white region', 'FontSize', fontSize);
ylabel('Gray Level', 'FontSize', fontSize);

[Kenni]

Many Thanks

[Image Analyst]

If that did the trick, please mark my answer as Accepted. Thanks in advance.

[Anders Mahler]

Hello, i have been looking at your exsample, you are writing about thresholding the EDT at 20, how would you do that? there is no threshold option in the bwdist() funktion?

[Image Analyst]

% Idenfity pixels farther away from any blob than 20
farAwayPixels = edtImage > 20; 
thresholdedEDT = edtImage; % Initialize
thresholdedEDT(farAwayPixels) = 0; % Mask it.