How to plot a 2D image profile along an arbitrary line passing through the image center

I would like to plot the profile of a 2D image alond lines passing from the image center. It is easier if the image is square (nxn). But it should be possible to do that for any image (nxm). For example the attached DICOM image is (768,1024). I read it through MatLab function 'dicomread'. I can plot profiles through the image center along the x-axis and the y-axis because I know the number of pixels along either direction. However, I would like to be able to plot a profile through the image center along the line joining the upper right and the lower left corner. In general I would like to plot a profile along a line, through the center, making an arbitrary angle with the x-axis. My problem is that I do not know how to calculate the number of pixels lying on arbitrarily slanted lines. Can you please help ?

Thank you very much.

Maura (<mailto:[e-mail] [e-mail]>)

Answer 1

[Image Analyst]

You can use improfile() and specify how many samples between the endpoints that you want. You can choose it so that the spacing is one pixel if you want. See my attached demo.

[Maura Monville]

You suggest to use 'improfile'. My problem is that I do not know how to find the endpoints. Given the arbitrary azimuthal angle 'phi' between the line and the x_axis I need to plot the image intensity along the intersection of such a line passing through the image center.

Such endpoints should be calculated as follows:

[xp yp] = size(img);
xmid = ceil(xp/2);
ymid = ceil(yp/2);
xi = xmid + radius * cos(phi);
yi = ymid + radius * sin(phi);
xf = xmid + radius * cos(phi + pi);
yf = ymid + radius * sin(phi + pi);

How far can I go along the line from the image center before crossing the image boundaries ? My images are not always square. That is, how long 'radius' can be ? If I can solve this problem then I will use MatLab command "improfile([xi yi], [xf yf])".

In your demo I cannot understand the following distances. Precisely the terms containing "1",

distanceToUL = sqrt((1-y)^2 + (1-x)^2)
distanceToUR = sqrt((1-y)^2 + (columns-x)^2)
distanceToLL = sqrt((rows-y)^2 + (1-x)^2)

Furthermore you seem to add the contribution of all pixels at the same distance from the chosen center (x,y). I do not have to estimate the average intensity along a circumference. I have to estimate the intensity along the chosen diameter. No averaging.

Can you please help ?

Thank you in advance.

maura

Maybe some interpolation is necessary, though.

[Image Analyst]

First of all, you made what is probably the #1 most common beginner mistake: thinking x,y is row,column. The size() function returns rows and then columns in that order. Thus you should not want to call rows xp since x implies the horizontal (columns) direction, and you should not want to call columns yp since y implies the vertical (rows) direction.

The distanceToUL, etc. terms are the distance from point (x,y) to the Upper Left, Upper Right, Lower Left, and Lower Right pixels. The (x,y) coordinates of those corners are (1,1), (columns, 1), (1, rows), and (columns, rows) respectively. The 1 is there because the index of the first row, or first column, is 1.

You should compute your end points as whatever they are, even if they're outside the image. If they're outside the image, you'll get NaN for the profile values. You can just ignore those. Look at this demo:

grayImage = imread('cameraman.tif');
subplot(2,1,1);
imshow(grayImage);
axis on;
x1 = -10
y1 = -40
x2 = 400
y2 = 511;
hold on;
plot([x1,x2], [y1,y2], 'b+-', 'LineWidth', 2);
[xi, yi, theProfile] = improfile(grayImage, [x1,x2], [y1,y2]);
subplot(2,1,2);
plot(theProfile, 'b*-');
grid on;
msgbox('Look at beginning and end of theProfile in the workspace');

Of course you'll need to keep track of what distances they are. If you don't specify the number of points then they're spaces a pixel apart. So index 1 is a distance of, say, 1 and index of 10 is a distance of 10. So if your non-NaN indexes go from 3 to 300, and indexes 1 and 2 and 301 and beyond were Nan, then your valid distances along the line go from 3 to 300. Of course you can shift those however you want depending on what pixel along the line you consider to be the origin. For example if element 103 was the center, then 3 would be a distance of -100, and 300 would be a distance of (300-103).

[Maura Monville]

Thank you. I attached the DICOM image I am using to test my script 'of10.mat'. I am pursuing the usage of improfile. I still have two problems: (1) finding the pixel indices of the two segments ends. I am not computing them right. In fact, if I merely plot the profile I get the picture 'Profile10_10.jpg' which is not simmetric with respect to the plot center. (2) I cannot superimpose the profile to the image. In fact if I use 'imshow' I get the following message: >> imshow(img) Warning: Image is too big to fit on screen; displaying at 67% > In images.internal.initSize (line 71) In imshow (line 305)

If I use 'imagesc' then I can display the image ('Of10_10.jpg') but plotting the profile after 'hold on' produced no profile at all ... I suspect there is a scaing factor difference between the picture and the profile .... ????

My code is in the attached file 'Radial_Profile.m'

Can you please help me ?

Thank you in advance, Maura

Answer 2

[John BG]

Maura

plot does not work if you give it theProfile directly. Either just give a single vector and plot will plot it, or give pairs of ref_vector_1 vector_1, ref_vector_2 vector_2, .. and plot will plot each vector along each respective reference.

To do so replace your line with last attempt plot with the following:

L11=theProfile(:,:,1)'
L12=theProfile(:,:,2)'
L13=theProfile(:,:,1)'
nL1=[1:1:length(L11)]
figure(2);plot(nL1,L11,'r','LineWidth',2)
figure(3);plot(nL1,L12,'g','LineWidth',2)
figure(4);plot(nL1,L13,'b','LineWidth',2)

If you directly combine Red Green Blue components you may lose grey scale details so critical in medical images.

I split the components so you see them, before building a Luminance Y approximating the Grey scale you sent me by email, which is what you probably want. Human eye gives more importance to green, and that translates into losing grey shades if not applying the right proportion. There is a standard called CCIR601 that already has the right coefficients to combine RGB without losing grey shades that would mean missing details in an X ray image.

profile_Y_CCIR601=.200*L11 +.587*L12  +.114*L13
figure(5);plot(profile_Y_CCIR601,'k','LineWidth',2),grid on

IMPORTANT:

Yet this is not the section of the image with medical resolution that you sent by email, just an approximation.

1.- RESOLUTION LOSS:

Please bear in mind that from the medical image you sent by email to the image you used to apply the Image Analyst script there is a loss of resolution because the original image had a wider pixel range than the uint8 of standard RGB [0 255].

Now when using the standard script supplied by the Image Analyst you show the same of what looks like an X ray aperture (confirmation would be appreciated), but the slot is magenta while the surrounding is yellow, and the pixels do not seem to belong to medical image range, because they are only within [0 255].

2.- DICOM ENCODING NOT CONSIDERED:

Just from a quick browsing of

http://dicom.nema.org/

questions for you:

Q1.- did you really use a DICOM viewer to send the aperture screen shot, or you just took a photo without much lighting, with a common off-the-shelf optic camera, or your mobile phone?

Q2.- do you use a DICOM viewer like the free one available here

http://dicom.nema.org/

I ask these 2 questions because the pixels of the initial image sent by email are encoded. May be you are already familiar with such encoding but understanding it means having gone through the standard or through a good book explaining such standard.

Again, without knowing DICOM encoding, approximating an X ray frame with an optic screen shot, seems quite a long shot to me.

Repeat using the original pixels to avoid losing resolution, but first get to know what each DICOM pixel really means.

If you find this answer of any help solving your question, please click on the thumbs-up vote link,

thanks in advance

John

[e-mail]

[Maura Monville]

I answered your questions to your persnal email address.

My problem is to be unable to get the image intensity profile along a diameter of the picture. I agree with you about the fusion of the 3 profiles. That is really interesting. Thank you. Nevertheless, have you notices that all the profiles you posted are not simmetrical ? This is because my code does not calculate correctly the pixel indices of the two diameter ends. Please, try out the attached script on the DICOM image I sent you (that here I have attached as a .mat file as .dcm is not allowed). Such a script generates the profiles along the x_axis and the y_axis. They are perfectly symmetrical with respect to the image center. How can I get a similar result plotting the profile along an arbitrary diameter of the picture ? Thanks. maura (<mailto:[e-mail] [e-mail]>)

[John BG]

i cropped the yellow magenta image (that is not a grey scale like i guess all DICOM are) to remove the indices along the axes.

Let me point out that the eccentricity of the previous cross section is not relevant, not even that the center point of all cross sections is really centered, at least not to solve the questions you have asked in the MATLAB forum that are mainly related to Image Processing, not to antenna beam forming, let me explain why I said this:

You are the X-ray expert, not me, but on one side you want to slice and image produced with X-ray illuminating on one direction, machine-to-patient.

But then you took a picture of the illuminating cavity?

Backwards patient-to-machine? with room lighting?

how else would you see the cavity aperture on the image you later supplied with a sample cross section.

It's my opinion that you should discard anything that is not the result of X-ray illumination, visible light and X-ray bands do not overlap, do they?

So, let's focus on correctly processing and getting the slicing of the imagery you want, produced for medical analysis.

Correct me if wrong, but the targets of your questioning are (amend as appropriate):

1.- build a volume to be able to use command slice

2.- use a rotating plane to slice around an axis parallel z axis that you want it right in the middle of the aperture of the X-ray slot antenna.

Please do the following:

1. rename the dicom_image.dcom to dicom_image_dcom.mat ignore any warning of your operative system, just replace the dot, and add '.mat' or '.m' at the tail of the file name.

2. attach the renamed original image to another comment.

3. With whatever DICOM viewer you use, try exporting the same image, the visual of the DICOM image you want to section, to another format like .jpeg or .bmp

I ask you for all this in order to see whether DICOM to RGB is really deteriorating that much the resolution that experienced doctors only achieve throughout years of X-ray radiography reading.

Regards

John [e-mail]