How can the major axis of a 3d object be determined?

I have some rodlike objects in 3d, below is an example. Is there a way to find its longitudinal axis? I'm trying to generate a bundle of such objects at different orientations and need to know the axis to be able to rotate them about that axis (and properly place them at the correct distance)

Answer 1

[Carson Purnell]

Seems like no matlab code has this functionality built in, so I ended up doing it myself. Demo code, assuming you already have a shape within 3d array vol:

%generate points from the volume
[X, Y, Z] = meshgrid(1:size(vol,2), 1:size(vol,1), 1:size(vol,3));
points = [X(:) Y(:) Z(:) vol(:)];
nonzero = points(:,4)>0; %threshold to retain points
filtered = points(nonzero,1:3);

%generate an orientation vector from points
cen = mean(filtered,1);
[~,~,V] = svd((filtered-cen),'econ'); %singular value decomposition, econ for speed
d = V(:,1)'; %the slope vector, ready for input to imrotate3

%demo of line
count =0;
for i=-50:.1:50 %show line
    count = count+1;
    lin(count,:) = d*i+cen;
end

close all
hold on
plot3(filtered(:,1),filtered(:,2),filtered(:,3),'.')
plot3(lin(:,1),lin(:,2),lin(:,3),'o');
axis equal

function is on the FEX, let me know if it works for anyone else.

https://www.mathworks.com/matlabcentral/fileexchange/113325-shapeaxis-compute-the-axis-orientation-of-3d-shapes

Answer 2

[Kevin Holly]

You can figure out the orientation with regionprops3 assuming you have a 3D matrix of the model.

Example:

sample = zeros(100,100,100);
sample(5,:,:)=eye(100);
SE = strel("disk", 3);
new_sample=imdilate(sample,SE);
isosurface(new_sample)
axis equal
xlim([0 100])
ylim([0 100])
zlim([0 100])

rp = regionprops3(new_sample,'Orientation');
rp.Orientation

ans = 1×3

-90.0000 45.0000 0

Let me rotate it with imrotate3

figure
new_sample_rotated = imrotate3(new_sample,90,rp.Orientation,'nearest','loose','FillValues',0);
isosurface(new_sample_rotated)
axis equal
xlim([0 100])
ylim([0 100])
zlim([0 100])

figure
new_sample_rotated = imrotate3(new_sample,35,rp.Orientation,'nearest','loose','FillValues',0);
isosurface(new_sample_rotated)
axis equal
xlim([0 100])
ylim([0 100])
zlim([0 100])

[Carson Purnell]

Thanks, that function looks promising. However, the orientation specification is not a vector description like I expected and have been using. An example: the original object (purple) is aligned almost perfectly to the Z axis, so any rotation of the object on that axis, [0 0 1], must create an overlapping object. The overlapping yellow shape shows that this works, but the separate yellow object is the one rotated along the regionprops3 output orienation, which is ~[55 90 0].

Now trying to figure out how to translate from that weird specification to an axis.

[Kevin Holly]

"Euler angles, returned as a 1-by-3 vector. The angles are based on the right-hand rule. regionprops3 interprets the angles by looking at the origin along the x-, y-, and z-axis representing roll, pitch, and yaw respectively. A positive angle represents a rotation in the counterclockwise direction. Rotation operations are not commutative so they must be applied in the correct order to have the intended effect."

sample = zeros(100,100,100);
sample(5,:,:)=eye(100);
SE = strel("disk", 3);
new_sample=imdilate(sample,SE);
isosurface(new_sample)
axis equal
xlim([0 100])
ylim([0 100])
zlim([0 100])

rp = regionprops3(new_sample,'Orientation','Centroid');
rp.Orientation

ans = 1×3

-90.0000 45.0000 0

rp.Centroid

ans = 1×3

50.5000 5.0000 50.5000

figure
% new_sample=imtranslate(new_sample,-rp.Centroid);
new_sample_rotated = imrotate3(new_sample,-45,[0 1 0]);
isosurface(new_sample_rotated)
axis equal
xlim([0 100])
ylim([0 100])
zlim([0 100])

rp = regionprops3(new_sample_rotated,'Orientation','Centroid');
rp.Orientation

ans = 1×3

90 0 0

rp.Centroid

ans = 1×3

71.5000 5.0000 71.5000

figure
% new_sample_rotated=imtranslate(new_sample,-rp.Centroid);
new_sample_rotated = imrotate3(new_sample_rotated,90,[0 0 1]);
isosurface(new_sample_rotated)
axis equal
xlim([0 100])
ylim([0 100])
zlim([0 100])

rp = regionprops3(new_sample_rotated,'Orientation');
rp.Orientation

ans = 1×3

0 0 0

Using rotate:

figure
isosurface(new_sample_rotated)
axis equal
xlim([0 100])
ylim([0 100])
zlim([0 100])
h=gca;
h.Children(2)

ans = 
  Patch with properties:

    FaceColor: 'flat'
    FaceAlpha: 1
    EdgeColor: 'none'
    LineStyle: '-'
        Faces: [5148×3 double]
     Vertices: [2576×3 double]

  Show all properties

rotate(h.Children(2),[1 0 0],45)

figure
isosurface(new_sample_rotated)
hold on
isosurface(new_sample_rotated)
axis equal
xlim([0 100])
ylim([0 100])
zlim([0 100])
h=gca;
rotate(h.Children(1),[1 0 0],35)

[Carson Purnell]

I've found no way to convert this transformation into the real axis of the object, even matlab's method in the aerospace toolbox does not generate anything close to the real axis.