Triangulating fiducial marker position in c-arm image (DICOM) using stereoParameters

Hi,

I have 2 DICOM images exported from a c-arm imaging device. The images are 60 degrees apart in rotation and there is no translation between the two images. There is a marker in each image of which I would like to return the 3-D coordinates in world coordinate system.

The source-imager distance, pixel size, coordinates of the marker in each image, and angle between the two images are all known.

I'm attempting to use stereoParameters and triangulate to achieve this however I am unable to change the FocalLength parameter manually.

Is anyone able to assist me with this? Or recommend an alternative method to achieve my aim?

Regards, Joel

[Matt J]

Yan Li commented:

Hi,Joel I am a graduate student now and interested in your C_arm image because I wanna study the parameters of them, could you upload them and share, I will appreciate it ~ THANKS ：） Leanne

Answer 1

[Matt J]

If you have the 3x4 projection matrices for the two gantry positions, then you can use the attached routine to do the triangulation. It isn't the optimal triangulation algorithm, but if your 2D localization errors are small, it serves most practical purposes.

[Joel Poder]

Hi Matt,

Thanks so much for the code. I'm new to this field and was hoping you could help further by assisting me in calculating my projection matrix.

My SID = 1650 mm, pixel size = 0.388 mm, centre pixel x = 512, centre pixel y = 384, source-origin distance = 1000 mm, angle between the two projections is 90 degreees, there is no translation.

I have calculated my intrinsic matrix to be

Intrinsic = [640.2, 0, 0; 0, 640.2, 0; 512, 384, 1]

Rotation matrix to be:

Rotation = [0, -1, 0; 1, 0, 0; 0, 0, 1]

Translation matrix to be:

Translation = [0; 1000; 0]

Are you able to please confirm this is correct format to be used as input into your function?

Thanks in advance, Joel

[Joel Poder]

Marker position coordinates in image 1

x,y = 524.7754, 372.8574

Marker position coordinates in image 2

x,y = 392.1953, 375.5078

The code is returning

X = [NaN; NaN; Inf]

I'm likely doing something wrong, but not sure what.

[Joel Poder]

Hi Matt,

I was hoping if you could please just clarify the units of the input data for the function. Should the origin coordinates in the intrinsic matrix be in pixel coordinates?

Should the source-origin distance in the translation matrix be in mm?

And finally, should the 2D projection coordinates (xx in your function) be pixel coordinates as well?

Thanks in advance...

Regards, Joel

[Matt J]

Hi Joel,

The coordinates are up to you, but they must be consistent. In other words all distances must be measured in the same units (either mm or pixels) in both the input xy coordinates and the intrinsic parameters.

Below, I illustrate the use of the code when all of the data you mentioned is converted to millimeters, but you can test for yourself that you get the same result when both the intrinsics and xy are re-expressed in pixels. I had to make some assumptions (see the A matrix) about how your 2D image axes are oriented relative to the 3D world coordinates.

pix=0.388; %pixel size

K = [1650, 0, 0; 0, 1650, 0; 512*pix, 384*pix, 1].' ;%intrinsic

C=[0;-1000;0]; %source location

A=[1 0 0 ; 0 0 -1; 0 1 0]; %assumed world coordinates of camera axes in position 1

R=[0, -1, 0; 1, 0, 0; 0, 0, 1]; %90 deg rotation

E1=A*[eye(3),-C ]; %Extrinsics position 1

E2=(A*R.')*[eye(3), -R*C]; %Extrinsics position 2

P(:,:,2)=K*E2; %projection matrix stack
P(:,:,1)=K*E1;

xy=pix*[524.7754, 372.8574; 392.1953, 375.5078].'; %observed 2D positions

xyz=triangulateRays(P,xy), %3D location in mm

This gives the result

xyz =

      2.9196
    -28.0961
      2.2691

as the triangulated 3D location.

[Joel Poder]

Hi Matt,

Thanks so much for all of this. I think I'm almost there. I've attached my code, where i've incorporated your triangulateRays.m

It is almost working however the code is failing on the line

Z=Pij-bsxfun(@times,xx,Pk)

The error is "Array dimensions must match for binary array op"

I've tried to fault find but cant seem to figure out why its not working.

Thanks in advance,

Joel

[Matt J]

I can't run your code without the input .dcm files, but I think this,

x(:,:,1) = [x1_mm, y1_mm];
x(:,:,2) = [x2_mm, y2_mm];

should really be this

x=[x1_mm, y1_mm; x2_mm, y2_mm].';

The fiducial (x,y) coordinates should be arranged in columns.

[Joel Poder]

Thanks Matt, that fixed the error. It's running now but I have one final question.

I'm identifying 3 fiducials on 2 different C-Arm images, however the code is returning a 3x6 matrix, but from my understanding it should be a 3x3 (i.e. M=3).

I've used the input for x as you suggested above and P is calculated as per the code I sent you.

Thanks for much for all of your help, its greatly appreciated.

Joel

[Joel Poder]

At the moment the input for x is a 2x6 Matrix, but I think it should be a 2x3x2 Matrix?

This is what I was attempting to do originally, but the error kept coming up in the code.

[Matt J]

Yes, it should be a 2x3x2 matrix, but xy(1,:,:) should be the x-coordinates and xy(2,:,:) should be the y-coordinates.

[Joel Poder]

Brilliant.

Thanks so much Matt