Finding coordinates of point(s) that are specific distance(s) away on a slanted straight line

Hello, I am a physician trying to locate points of interest on medical images. My problem boils down to the following:

I have an image with 7 points of interest that I need to identify on my image. The 7 points lie in a slanted straight line like so:

I could take advantage of image properties to find the coordinates of point A and C. d1, d2, and d3 are all known values. What script should I run to find the coordinates of the 2 crosses d1 distance from A? And similarly, how do I find the coordinates for B, and the 2 coordinates d3 distance from that?

Thank you so much for your help!

[Star Strider]

Do you have any recommendations for radiology text or reference books relevant to both clinical and technical perspectives? (I’d appreciate Author-Title-ISBN so I know I’m getting the correct ones.) I’ve been looking for a while, but haven’t found any that seem to cover what I want. You’re the perfect person to ask!

Answer 1

[Star Strider]

Not a function, but a short bit of code:

A = [5 10];                                                 % Create Data
C = [7 15];                                                 % Create Data
d1 = 2.1;
d2 = 3.8;
d3 = 1.5;

x = [A(1); C(1)];                                           % Define Independent Variables
y = [A(2); C(2)];                                           % Define Dependent Variables

P = [ones(size(x(:)))  x(:)]\y(:);                          % Calculate Intercept, Slope For Line

phi = atan2(diff(y), diff(x));                              % Slope Angle (rad)

LineDist = @(V,d,phi) [V(1)+d*cos(phi)  V(2)+d*sin(phi);  V(1)-d*cos(phi)  V(2)-d*sin(phi)];

d1X = LineDist(A,d1,phi);                                   % ‘d1’ Coordinates
B   = LineDist(C,d2,phi);                                   % ‘B’ Coordinates
d3X = LineDist(B(1,:),d3,phi);                              % ‘d3’ Coordinates

xv = linspace(min(x)-3, max(x)+3);                          % X-Vector For Plot
yv = [ones(size(xv(:)))  xv(:)]*P;                          % Y-Vector For Plot

figure(1)
plot(xv, yv)
hold on
plot(x, y, 'ok', 'MarkerFaceColor','k')
plot(d1X(1,1), d1X(1,2), 'xr')
plot(d1X(2,1), d1X(2,2), 'xr')
plot(B(1,1), B(1,2), 'xk')
plot(d3X(1,1), d3X(1,2), 'xr')
plot(d3X(2,1), d3X(2,2), 'xr')
hold off
grid
text(A(1)+0.1,A(2), 'A', 'HorizontalAlignment','left')
text(B(1,1)+0.1,B(1,2), 'B', 'HorizontalAlignment','left')
text(C(1)+0.1,C(2), 'C', 'HorizontalAlignment','left')
text(d1X(:,1)+0.1,d1X(:,2), 'd1', 'HorizontalAlignment','left')
text(d3X(:,1)+0.1,d3X(:,2), 'd3', 'HorizontalAlignment','left')

The plot:

Board-Certified Internist here, with M.Sc. in Biomedical Engineering!

[Brian]

Thank you so much for your help. Let me try to digest this. I have another quick question:

I have the slopes (calculated by y1-y0/x1-x0) already. How could I incorporate that information into your code? In other words, how do I quickly convert that to phi?

Radiologist in the training here! Sadly I have no prior engineering background and is just picking up Matlab.

[Star Strider]

Thank you! My pleasure!

I hope postgraduate clinical training is more humane now than it was when I went through it. I learned FORTRAN as an undergraduate Chemistry major, then after getting my Boards, went back to graduate school to understand instrumentation, signal processing, and control theory. Picked up MATLAB along the way in further postgraduate training. Retired now, but keep up with IM and BME.

If you want to do image processing with MATLAB, ‘Image Analyst’ is our local expert here on MATLAB Answers, so feel free to post any questions about image processing (preferably with the Image Processing Toolbox) you may have. (He also has a collection of demos in the File Exchange.) It always helps to upload the image you want help with. Use the brown-framed green ‘picture’ icon. It will appear the same way my plot does here (saved as a .png image first).

I’ll help as I can.

[Brian]

Clinical training is by the hose per usual, but likely more manageable now than before with more accessible information and work hour regulations. Exciting to hear that you did a rather unique track! Very cool!

I dabble with image processing projects every so often and Image Analyst did help me tremendously in the past. This is a wonderful community indeed.

My silly question still stands: I have slope (derived from A, C in rise/run format). How could I convert that to the phi that you came up with? Asking so I could use the very helpful magical line of:

LineDist = @(V,d,phi) [V(1)+d*cos(phi) V(2)+d*sin(phi); V(1)-d*cos(phi) V(2)-d*sin(phi)];

I feel like I have all the values in place, just don't know how to put them together/convert/apply to get my final results.

[Star Strider]

I envy your work hour regulations. Our 100-hour weeks were the rule rather than the exception back then. The best day of my residency was the last one! Going back and getting my M.Sc was a way to restore my sanity.

I wrote my code so you could just input the coordinates of ‘A’, ‘C’ and the three distances, and everything else will take care of itself. It does all the calculations I anticipate your needing. My ‘phi’ is simply the quadrant-corrected arctangent (using atan2) of the rise/run. If you already have the slope in rise/run form, ‘phi’ is the arctangent (using atan) of the slope you calculated. I kept all the angles in radian rather than degree units. (The slope in the more traditional sense is ‘P(2)’, and then ‘phi=atan(P(2))’.)

My ‘LineDist’ anonymous function (see the section on them in Function Basics for details if you are not familiar with them) takes any point ‘V’ defined as the (1x2) vector of coordinates ‘[Vx Vy]’, distance ‘d’ and regression-line-slope-angle ‘phi’ and uses fundamental trigonometry to calculate the coordinates of the points ‘d’ distance on both sides of ‘V’ on the line.

I could easily write my code as a function file (without the plot) if that would make it easier for you to use. It would take the (1x2) vectors ‘A’, ‘C’ and the three scalar distances as inputs, and return ‘d1X’, ‘B(1,:)’, and ‘d3X’ as outputs.

[Brian]

You must be well aware of the 80-hour per week regulations now. Though I think with the increased demand of documentations, the 80 hours doesn't quite seem enough. I am glad you enjoyed your M. Sc!

I didn't see you comments in time but did some reading to make atan2 work beautifully along with your LineDist function. Thanks so much for your very detailed explanations.

I actually do have another point (say E), that is lies on a line perpendicular to AC, passes by C, at given distance d4 from C. I am still trying to figure out how to take advantage of the atan2 results of AC to find the radian for E (I can hear my middle school geometry teacher cry in the background...but math is such a distant memory).

And just for kicks, I am using the script to get Matlab to automate CT ROI measurements to analyze our instruments :)

[Star Strider]

I’ve read about them, but 80 hours is still a lot. At least it should give you the opportunity to read. We tried to keep up, but as soon as we started reading anything, we’d go immediately to sleep. I finally caught up with both sleep and reading in my Endocrine fellowship. I then chose to pursue engineering rather than take my Endocrine Boards, since my IM Boards were enough. Engineering school was actually more difficult than Med school, but it was also more fun because the feedback was immediate. The IM Boards themselves were the most difficult exams I ever sat.

If I remember correctly, a line orthogonal to a given line (in a plane) has a slope (not the angle of the slope) that is the negative of that of the given line. (I proved that while back using the dot-product, but I’m not yet awake enough to prove it again this morning.) So, with one point (apparently ‘C’) and a known slope of the original line, calculating the equation of the orthogonal line is straightforward.

I didn’t test this, but you could also use ‘LineDist’ with the coordinates of ‘C’ (for ‘V’), the distance, and phi±(pi/2) for ‘phi’, and then choose the appropriate row (the coordinates it returns should be symmetric on each side of the line).

It didn’t occur to me when I wrote it that ‘LineDist’ could calculate the orthogonal line as well, but it should be able to. Give that a go first.

[Image Analyst]

"a line orthogonal to a given line (in a plane) has a slope (not the angle of the slope) that is the negative of that of the given line" <== it's slope is actually -1/slope, the negative inverse.

[Star Strider]

The ‘phi±(pi/2)’ idea works. It’s necessary to use axis equal on the plot to see it correctly.

Adding the off-line distance ‘d4’ and calculating the perpendicular line at ‘C’ as the ‘PrpC’ coordinate matrix (plotted in green here), the complete code becomes:

A = [5 10];                                                 % Create Data
C = [7 15];                                                 % Create Data
d1 = 2.1;
d2 = 3.8;
d3 = 1.5;

x = [A(1); C(1)];                                           % Define Independent Variables
y = [A(2); C(2)];                                           % Define Dependent Variables

P = [ones(size(x(:)))  x(:)]\y(:);                          % Calculate Intercept, Slope For Line

phi = atan2(diff(y), diff(x));                              % Slope Angle (rad)

LineDist = @(V,d,phi) [V(1)+d*cos(phi)  V(2)+d*sin(phi);  V(1)-d*cos(phi)  V(2)-d*sin(phi)];

d1X = LineDist(A,d1,phi);                                   % ‘d1’ Coordinates
B   = LineDist(C,d2,phi);                                   % ‘B’ Coordinates
d3X = LineDist(B(1,:),d3,phi);                              % ‘d3’ Coordinates

d4 = 5;
PrpC = LineDist(C,d4,phi+pi/2);

xv = linspace(min(x)-3, max(x)+3);                          % X-Vector For Plot
yv = [ones(size(xv(:)))  xv(:)]*P;                          % Y-Vector For Plot

figure(1)
plot(xv, yv)
hold on
plot(x, y, 'ok', 'MarkerFaceColor','k')
plot(d1X(1,1), d1X(1,2), 'xr')
plot(d1X(2,1), d1X(2,2), 'xr')
plot(B(1,1), B(1,2), 'xk')
plot(d3X(1,1), d3X(1,2), 'xr')
plot(d3X(2,1), d3X(2,2), 'xr')
plot(PrpC(:,1)', PrpC(:,2)', 'p-g')
hold off
grid
text(A(1)+0.1,A(2), 'A', 'HorizontalAlignment','left')
text(B(1,1)+0.1,B(1,2), 'B', 'HorizontalAlignment','left')
text(C(1)+0.1,C(2), 'C', 'HorizontalAlignment','left')
text(d1X(:,1)+0.1,d1X(:,2), 'd1', 'HorizontalAlignment','left')
text(d3X(:,1)+0.1,d3X(:,2), 'd3', 'HorizontalAlignment','left')
text(PrpC(:,1)+0.1,PrpC(:,2), 'd4', 'VerticalAlignment','bottom')
axis equal

with the plot:

Answer 2

[Image Analyst]

Knowing points A and C, the equation of the line is

slope = (yc-ya) / (xc - xa);
y = slope * (x - xa) + ya;

To get the x for the left and right d1 point from A

angle = atan(slope);
xd1Left = xa - d1 * cos(angle);
xd1Right = xa + d1 * cos(angle);

To get point B from known point C

xb = xc + d2 * cos(angle);

To get the d3 points

xd3Left = xc + (d2-d3) * cos(angle);
xd3Right = xc + (d2+d3) * cos(angle);

To get the y values for any of those x values, plug them in for x in the formula for the line:

y = slope * (x - xa) + ya;