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

11 views (last 30 days)
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)

Accepted Answer

Carson Purnell
Carson Purnell on 15 Jun 2022
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.

More Answers (1)

Kevin Holly
Kevin Holly on 9 Jun 2022
Edited: Kevin Holly on 9 Jun 2022
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])
  3 Comments
Kevin Holly
Kevin Holly on 10 Jun 2022
Edited: Kevin Holly on 10 Jun 2022
"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
Carson Purnell on 13 Jun 2022
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.

Sign in to comment.

Categories

Find more on 2-D and 3-D Plots in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!