How can I retrieve x,y coordinates of a sliced plane through a 3D volume using the obliqueslice function?

Hello,

I am using the function obliqueslice to get an aribtrary (not orthogonal) sliced plane of my 3D dataset. I want to get the coordinates of the point that used as input in the resulting image slice. It is necessary to know how the slice is build up. The documentation of the function inlcudes a brief description of the orientation of slice in image plane, but I don't understand exactly how the output image is build up.

My question is: what is meant by "The slice corner that lies close to the origin of the last slice (0,0,P), in the volumetric data constitutes the upper-left pixel in the image plane."?

My interpretation of this explanation is that the slice corner with the shortest Euclidean distance to the point x=0,y=0 and z=size(Volume,3) is always in the left upper corner of the 2D image. However, this is not true in all cases. For example (with example code from documentation):

load mri
V = squeeze(D);

point = [73 5 15.5];
normal = [5 5 -20];
[B,x,y,z] = obliqueslice(V,point,normal);

fig = figure;
axes1 = axes('Parent',fig,'Position',[0.13 0.4 0.37 0.5]);
surf(x,y,z,B,'EdgeColor','None','HandleVisibility','off','Parent',axes1);
grid on
view([-62.6 13.8])
colormap(gray)
xlabel('x-axis')
ylabel('y-axis');
zlabel('z-axis');
title('Position of Slice in 3-D Coordinate Space')

hold on
plot3(point(1),point(2),point(3),'or','MarkerFaceColor','r');
plot3(normal(1),normal(2),normal(3),'ob','MarkerFaceColor','b');
plot3(0,0,size(V,3),'oy','MarkerFaceColor','y');
plot3(x(1,1),y(1,1),z(1,1),'om','MarkerFaceColor','m');
plot3(x(1,end),y(1,end),z(1,end),'oc','MarkerFaceColor','c');
plot3(x(end,end),y(end,end),z(end,end),'og','MarkerFaceColor','g');
plot3(x(end,1),y(end,1),z(end,1),'ok','MarkerFaceColor','k');
hold off
legend('Point in the volume','Normal vector','Origin last slice','1','2','3','4','Position',[0.15 0.1 0.3 0.08])

axes2 = axes('Parent',fig,'Position',[0.59 0.39 0.33 0.6]);
imshow(B,[],'Parent',axes2)
title('Output Slice in Image Plane')

corner1 = [x(1,1),y(1,1),z(1,1)];
corner2 = [x(1,end),y(1,end),z(1,end)];
corner3 = [x(end,end),y(end,end),z(end,end)];
corner4 = [x(end,1),y(end,1),z(end,1)];
o_lastslice = [0, 0, size(V,3)];

length_corner1 = sqrt((corner1(1)-o_lastslice(1))^2+(corner1(2)-o_lastslice(2))^2+...
(corner1(3)-o_lastslice(3))^2);
length_corner2 = sqrt((corner2(1)-o_lastslice(1))^2+(corner2(2)-o_lastslice(2))^2+...
(corner2(3)-o_lastslice(3))^2);
length_corner3 = sqrt((corner3(1)-o_lastslice(1))^2+(corner3(2)-o_lastslice(2))^2+...
(corner3(3)-o_lastslice(3))^2);
length_corner4 = sqrt((corner4(1)-o_lastslice(1))^2+(corner4(2)-o_lastslice(2))^2+...
(corner4(3)-o_lastslice(3))^2);

Here, corner3 has the shortest distance to the origin of the last slice, however, corner 1 is taken as upper-left pixel in the image plane.

I am looking forward to get some more information about this function obliqueslice. Other ideas how to get the coordinates of the point in the 2D slice are also welcome. Many thanks!

Answer 1

[Boubat Matthieu]

Hello,

I don't know anything about the upper-left pixel of the slice, but this is how I find a point's indices in a slice:

load mri
V = squeeze(D);

pointSubscriptIndicesInVolume = [40 60 13];
normal = [5 5 -20];
[B,x,y,z] = obliqueslice(V,pointSubscriptIndicesInVolume,normal);

pointLinearIndexInSlice = find( round(x) == pointSubscriptIndicesInVolume(1) &...
                                round(y) == pointSubscriptIndicesInVolume(2) &...
                                round(z) == pointSubscriptIndicesInVolume(3) );
                      
[pointRow,pointColumn] = ind2sub(size(B),pointLinearIndexInSlice);

This solution requires to use the point's subscript indices instead of its coordinates.

I haven't checked that it will always provide a solution, and only one solution, but it has not failed me so far.

[Astrid]

Hello Boubat Matthieu,

Many thanks for your alternative solution. It works!

I din't always obtain a result in my data, so I'm now using ceil() instead of round() and I'm not only looking at the the pixel equal to the point's subscript indices, but also to one next to it. This sometimes give more than one solution, so I pick the first one (random, but in my case always close enough to the original coordinates):

pointLinearIndexInSlice = find(ceil(x)==ceil(pointSubscriptIndicesInVolume(1)) &...
                               ceil(y)==ceil(pointSubscriptIndicesInVolume(2)) &...
                               ceil(z)==ceil(pointSubscriptIndicesInVolume(3)) );
if isempty(pointLinearIndexInSlice)
    pointLinearIndexInSlice = find((ceil(x)==ceil(pointSubscriptIndicesInVolume(1))|ceil(x)==(ceil(pointSubscriptIndicesInVolume(1))-1)) &...
                                   (ceil(y)==ceil(pointSubscriptIndicesInVolume(2))|ceil(y)==(ceil(pointSubscriptIndicesInVolume(2))-1)) &...
                                   (ceil(z)==ceil(pointSubscriptIndicesInVolume(3))|ceil(z)==(ceil(pointSubscriptIndicesInVolume(3))-1)));   
end
    
for j=1:size(pointLinearIndexInSlice,1)
    [pointRow,pointColumn] = ind2sub(size(B),pointLinearIndexInSlice);
        if size(pointRow,2) > 1
            pointRow = pointRow(1);
            pointColumn = pointColumn(1);
        end
end