Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Fitting a curve into a 3D point cloud

Subject: Fitting a curve into a 3D point cloud

From: Stan

Date: 29 Sep, 2010 14:58:05

Message: 1 of 10

I have data that is made up of a 3d point cloud (x,y,z coordinates). It essentially looks like a curved tube. My problem is that I can't find a way to calculate the center-line of this "tube of points". I can only find information on fitting surfaces to 3d point clouds using least square methods.

Any ideas?

Help would be greatly appreciated!

Subject: Fitting a curve into a 3D point cloud

From: Matt J

Date: 29 Sep, 2010 16:00:18

Message: 2 of 10

"Stan " <zufallsacc@yahoo.com> wrote in message <i7vk5t$pir$1@fred.mathworks.com>...
> I have data that is made up of a 3d point cloud (x,y,z coordinates). It essentially looks like a curved tube. My problem is that I can't find a way to calculate the center-line of this "tube of points". I can only find information on fitting surfaces to 3d point clouds using least square methods.
>



%Simulated data on the parametric line [1;1;1]+t*[0 0 1]

N=10;
XYZ=bsxfun(@plus, [ 1 1 1], (-N:N).'*[0 0 1]);
 XYZ=randn(size(XYZ))*.01 + XYZ; %Add noise

 
%Fit to a line xyz0 + t*direction


xyz0=mean(XYZ);
A=bsxfun(@minus,XYZ,xyz0); %center the data


[U,S,V]=svd(A);

xyz0,
direction=cross(V(:,end),V(:,end-1) ),

Subject: Fitting a curve into a 3D point cloud

From: Matt J

Date: 29 Sep, 2010 16:59:06

Message: 3 of 10

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i7vnqi$rus$1@fred.mathworks.com>...

> %Fit to a line xyz0 + t*direction
>
>
> xyz0=mean(XYZ);
> A=bsxfun(@minus,XYZ,xyz0); %center the data
>
>
> [U,S,V]=svd(A);
>
> xyz0,
> direction=cross(V(:,end),V(:,end-1) ),
=================


Forget it. I hadn't noticed in your post that the tube was curved. Do you have a model for the shape of the tube, i.e., an equation it is supposed to obey?

Subject: Fitting a curve into a 3D point cloud

From: Stan

Date: 30 Sep, 2010 06:05:20

Message: 4 of 10

"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i7vr8q$f0m$1@fred.mathworks.com>...
> "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i7vnqi$rus$1@fred.mathworks.com>...
>
> > %Fit to a line xyz0 + t*direction
> >
> >
> > xyz0=mean(XYZ);
> > A=bsxfun(@minus,XYZ,xyz0); %center the data
> >
> >
> > [U,S,V]=svd(A);
> >
> > xyz0,
> > direction=cross(V(:,end),V(:,end-1) ),
> =================
>
>
> Forget it. I hadn't noticed in your post that the tube was curved. Do you have a model for the shape of the tube, i.e., an equation it is supposed to obey?

Thanks for your reply.

No, I don't have a model right now. I thought about just projecting the points onto the xy, xz and yz planes and then do 3 2D curve fits with polynomials of a sufficiently high order and thus get an approximation for the center line shape, but I guess this is not really elegant..and won't work in more complicated cases involving loops..

Subject: Fitting a curve into a 3D point cloud

From: John D'Errico

Date: 30 Sep, 2010 09:17:08

Message: 5 of 10

"Stan " <zufallsacc@yahoo.com> wrote in message <i819b0$rg6$1@fred.mathworks.com>...
> "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i7vr8q$f0m$1@fred.mathworks.com>...
> > "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i7vnqi$rus$1@fred.mathworks.com>...
> >
> > > %Fit to a line xyz0 + t*direction
> > >
> > >
> > > xyz0=mean(XYZ);
> > > A=bsxfun(@minus,XYZ,xyz0); %center the data
> > >
> > >
> > > [U,S,V]=svd(A);
> > >
> > > xyz0,
> > > direction=cross(V(:,end),V(:,end-1) ),
> > =================
> >
> >
> > Forget it. I hadn't noticed in your post that the tube was curved. Do you have a model for the shape of the tube, i.e., an equation it is supposed to obey?
>
> Thanks for your reply.
>
> No, I don't have a model right now. I thought about just projecting the points onto the xy, xz and yz planes and then do 3 2D curve fits with polynomials of a sufficiently high order and thus get an approximation for the center line shape, but I guess this is not really elegant..and won't work in more complicated cases involving loops..

No. That will work terribly. And high order polynomials
are not a good choice for the fit.

Are you willing to assume the tube has constant radius,
or can that vary? Is the tube perfectly circular in cross-
section at any point, if you look orthogonally to the
tube axis?

These are things you need to specify, if you will choose
a method for deriving a "model" for this.

John

Subject: Fitting a curve into a 3D point cloud

From: Stan

Date: 30 Sep, 2010 09:41:03

Message: 6 of 10

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <i81kik$mkp$1@fred.mathworks.com>...
> "Stan " <zufallsacc@yahoo.com> wrote in message <i819b0$rg6$1@fred.mathworks.com>...
> > "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i7vr8q$f0m$1@fred.mathworks.com>...
> > > "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i7vnqi$rus$1@fred.mathworks.com>...
> > >
> > > > %Fit to a line xyz0 + t*direction
> > > >
> > > >
> > > > xyz0=mean(XYZ);
> > > > A=bsxfun(@minus,XYZ,xyz0); %center the data
> > > >
> > > >
> > > > [U,S,V]=svd(A);
> > > >
> > > > xyz0,
> > > > direction=cross(V(:,end),V(:,end-1) ),
> > > =================
> > >
> > >
> > > Forget it. I hadn't noticed in your post that the tube was curved. Do you have a model for the shape of the tube, i.e., an equation it is supposed to obey?
> >
> > Thanks for your reply.
> >
> > No, I don't have a model right now. I thought about just projecting the points onto the xy, xz and yz planes and then do 3 2D curve fits with polynomials of a sufficiently high order and thus get an approximation for the center line shape, but I guess this is not really elegant..and won't work in more complicated cases involving loops..
>
> No. That will work terribly. And high order polynomials
> are not a good choice for the fit.
>
> Are you willing to assume the tube has constant radius,
> or can that vary? Is the tube perfectly circular in cross-
> section at any point, if you look orthogonally to the
> tube axis?
>
> These are things you need to specify, if you will choose
> a method for deriving a "model" for this.
>
> John

Thanks for the input, I was already suspecting that it might be a terrible idea ;).

Yes the tube has a given radius and is circular in cross-section.

Another idea I have is partitioning the tube and using RCA on each section, but I wonder if there is a better solution

Stan

Subject: Fitting a curve into a 3D point cloud

From: Sven

Date: 12 Oct, 2010 23:43:04

Message: 7 of 10

Hi all,

I have basically the same problem, perhaps slightly simpler in that:

1. My pts are of a curve with noise in 3D (rather than surface pts of a cylinder, but it's a similar idea).
2. I have known anchor pts at either end (ie, the ends of my curve are constrained)
3. My curve is of reasonably known shape (human rib - basically a skewed "C")
4. I only want a "reasonable" approximation.

For my reasonable approximation, I'm thinking of the following rules:
1. The final curve should be of, say, 10 equally spaced control points
2. The maximum angle between adjacent curve control points is, say, 30 degrees

I was thinking of approaching the problem by starting with 10 equally spaced control points from my start-to-end anchor pts. Then I would:

1. Assign each of my sample pts to its nearest control pt
2. Shift each control pt to the mean of the sample pts in its domain. Any control pts without sample pts in its domain get linearly distributed between its neighbours which do.
3. Re-sample the control points to be equally spaced
4. Repeat from step 1 until a certain number of iterations are reached

Any thoughts on this? Did the OP manage to get towards a solution for their problem?

Thanks,
Sven.

Subject: Fitting a curve into a 3D point cloud

From: Bruno Luong

Date: 13 Oct, 2010 00:36:03

Message: 8 of 10

"Sven" <sven.holcombe@gmail.deleteme.com> wrote in message <i92rq8$80r$1@fred.mathworks.com>...
> Hi all,
>
> I have basically the same problem, perhaps slightly simpler in that:
>
> 1. My pts are of a curve with noise in 3D (rather than surface pts of a cylinder, but it's a similar idea).
> 2. I have known anchor pts at either end (ie, the ends of my curve are constrained)
> 3. My curve is of reasonably known shape (human rib - basically a skewed "C")
> 4. I only want a "reasonable" approximation.
>
> For my reasonable approximation, I'm thinking of the following rules:
> 1. The final curve should be of, say, 10 equally spaced control points
> 2. The maximum angle between adjacent curve control points is, say, 30 degrees
>
> I was thinking of approaching the problem by starting with 10 equally spaced control points from my start-to-end anchor pts. Then I would:
>
> 1. Assign each of my sample pts to its nearest control pt
> 2. Shift each control pt to the mean of the sample pts in its domain. Any control pts without sample pts in its domain get linearly distributed between its neighbours which do.
> 3. Re-sample the control points to be equally spaced
> 4. Repeat from step 1 until a certain number of iterations are reached
>
> Any thoughts on this? Did the OP manage to get towards a solution for their problem?
>
> Thanks,
> Sven.

Why not simply fitting the curve with three parametrized splines, each applied respectively on x, y, and z coordinates?

Bruno

Subject: Fitting a curve into a 3D point cloud

From: Sven

Date: 13 Oct, 2010 15:08:03

Message: 9 of 10

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <i92utj$q9e$1@fred.mathworks.com>...
> "Sven" <sven.holcombe@gmail.deleteme.com> wrote in message <i92rq8$80r$1@fred.mathworks.com>...
> > Hi all,
> >
> > I have basically the same problem, perhaps slightly simpler in that:
> >
> > 1. My pts are of a curve with noise in 3D (rather than surface pts of a cylinder, but it's a similar idea).
> > 2. I have known anchor pts at either end (ie, the ends of my curve are constrained)
> > 3. My curve is of reasonably known shape (human rib - basically a skewed "C")
> > 4. I only want a "reasonable" approximation.
> >
> > For my reasonable approximation, I'm thinking of the following rules:
> > 1. The final curve should be of, say, 10 equally spaced control points
> > 2. The maximum angle between adjacent curve control points is, say, 30 degrees
> >
> > I was thinking of approaching the problem by starting with 10 equally spaced control points from my start-to-end anchor pts. Then I would:
> >
> > 1. Assign each of my sample pts to its nearest control pt
> > 2. Shift each control pt to the mean of the sample pts in its domain. Any control pts without sample pts in its domain get linearly distributed between its neighbours which do.
> > 3. Re-sample the control points to be equally spaced
> > 4. Repeat from step 1 until a certain number of iterations are reached
> >
> > Any thoughts on this? Did the OP manage to get towards a solution for their problem?
> >
> > Thanks,
> > Sven.
>
> Why not simply fitting the curve with three parametrized splines, each applied respectively on x, y, and z coordinates?
>
> Bruno

I don't think that I can do this immediately on my input points as aren't necessarily ordered and since the curve can be "C" shaped, they can double back in any given direction, so ordering by their coordinates can give unpredictable results. Below is some sample input - if I'm wrong and it's actually quite simple to fit 3 splines, can you show me how? Perhaps a good way to order my kind of input points would be by their distance from the start/end locations. I'm working on a function of the form
function cPts = fitCurveTo3DPts(XYZ, stPt, endPt, n)

I'll post here if/when it's done, but I'd appreciate any input you have.

Thanks,
Sven.

stPt = [53 136 -222];
endPt = [150 265 -112];
n = 10;
XYZ = [122.6333 277.2984 -113.9278
  125.6526 275.7797 -113.9278
  127.1622 275.7797 -113.9278
  127.1622 275.7797 -111.9149
  128.6719 274.2609 -111.9149
  130.1815 274.2609 -111.9149
  134.7105 271.2234 -111.9149
  137.7298 271.2234 -111.9149
  139.2394 271.2234 -111.9149
  149.8070 265.1484 -111.9149
  143.7684 268.1859 -109.9021
   44.1314 149.7234 -214.5722
   41.1121 154.2797 -210.5464
   39.6025 155.7984 -208.5335
   41.1121 154.2797 -208.5335
   38.0928 160.3547 -206.5206
   36.5832 172.5047 -198.4691
   39.6025 184.6547 -186.3918
   41.1121 189.2109 -182.3660
   41.1121 190.7297 -180.3531
   42.6218 193.7672 -178.3402
   42.6218 196.8047 -176.3273
   42.6218 199.8422 -174.3144
   44.1314 202.8797 -172.3015
   44.1314 204.3984 -170.2887
   45.6411 205.9172 -170.2887
   45.6411 208.9547 -168.2758
   45.6411 210.4734 -166.2629
   92.4403 271.2234 -128.0180
   95.4596 272.7422 -126.0052
   96.9692 272.7422 -123.9923
   98.4789 274.2609 -123.9923
   99.9885 275.7797 -123.9923
  101.4982 275.7797 -121.9794
  104.5175 275.7797 -121.9794
  107.5368 277.2984 -121.9794
  107.5368 275.7797 -119.9665
  109.0464 277.2984 -119.9665
  110.5561 277.2984 -117.9536
  112.0657 278.8172 -117.9536
  113.5754 280.3359 -115.9407
   69.7955 254.5172 -138.0825
   72.8148 257.5547 -138.0825
   71.3051 257.5547 -136.0696
   72.8148 259.0734 -136.0696
   75.8341 260.5922 -136.0696
   75.8341 260.5922 -134.0567
   78.8534 262.1109 -134.0567
   80.3631 263.6297 -134.0567
   78.8534 263.6297 -132.0438
   80.3631 265.1484 -132.0438
   83.3824 266.6672 -132.0438
   83.3824 266.6672 -130.0309
   84.8920 268.1859 -130.0309
   87.9113 269.7047 -130.0309
   86.4017 268.1859 -128.0180
   87.9113 269.7047 -128.0180
   90.9306 271.2234 -128.0180
   92.4403 271.2234 -126.0052
   54.6990 231.7359 -152.1727
   56.2086 233.2547 -150.1598
   57.7183 236.2922 -150.1598
   57.7183 237.8109 -148.1469
   59.2279 239.3297 -146.1340
   60.7376 240.8484 -146.1340
   60.7376 242.3672 -144.1211
   62.2472 243.8859 -144.1211
   62.2472 245.4047 -142.1082
   65.2665 248.4422 -142.1082
   66.7762 249.9609 -140.0954
   68.2858 252.9984 -140.0954
   68.2858 252.9984 -138.0825
   47.1507 213.5109 -164.2500
   48.6604 216.5484 -162.2371
   48.6604 218.0672 -160.2242
   50.1700 219.5859 -160.2242
   50.1700 221.1047 -158.2113
   51.6797 222.6234 -158.2113
   51.6797 224.1422 -156.1985
   53.1893 225.6609 -156.1985
   53.1893 228.6984 -154.1856
   38.0928 172.5047 -198.4691
   39.6025 172.5047 -196.4562
   39.6025 174.0234 -194.4433
   50.1700 139.0922 -224.6366
   53.1893 136.0547 -224.6366
   50.1700 139.0922 -222.6237
   53.1893 136.0547 -222.6237
   50.1700 139.0922 -220.6108
   51.6797 139.0922 -218.5979];

Subject: Fitting a curve into a 3D point cloud

From: Sven

Date: 13 Oct, 2010 15:56:04

Message: 10 of 10

"Sven" <sven.holcombe@gmail.deleteme.com> wrote in message <i94i0j$beo$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <i92utj$q9e$1@fred.mathworks.com>...
> > "Sven" <sven.holcombe@gmail.deleteme.com> wrote in message <i92rq8$80r$1@fred.mathworks.com>...
> > > Hi all,
> > >
> > > I have basically the same problem, perhaps slightly simpler in that:
> > >
> > > 1. My pts are of a curve with noise in 3D (rather than surface pts of a cylinder, but it's a similar idea).
> > > 2. I have known anchor pts at either end (ie, the ends of my curve are constrained)
> > > 3. My curve is of reasonably known shape (human rib - basically a skewed "C")
> > > 4. I only want a "reasonable" approximation.
> > >
> > > For my reasonable approximation, I'm thinking of the following rules:
> > > 1. The final curve should be of, say, 10 equally spaced control points
> > > 2. The maximum angle between adjacent curve control points is, say, 30 degrees
> > >
> > > I was thinking of approaching the problem by starting with 10 equally spaced control points from my start-to-end anchor pts. Then I would:
> > >
> > > 1. Assign each of my sample pts to its nearest control pt
> > > 2. Shift each control pt to the mean of the sample pts in its domain. Any control pts without sample pts in its domain get linearly distributed between its neighbours which do.
> > > 3. Re-sample the control points to be equally spaced
> > > 4. Repeat from step 1 until a certain number of iterations are reached
> > >
> > > Any thoughts on this? Did the OP manage to get towards a solution for their problem?
> > >
> > > Thanks,
> > > Sven.
> >
> > Why not simply fitting the curve with three parametrized splines, each applied respectively on x, y, and z coordinates?
> >
> > Bruno
>
> I don't think that I can do this immediately on my input points as aren't necessarily ordered and since the curve can be "C" shaped, they can double back in any given direction, so ordering by their coordinates can give unpredictable results. Below is some sample input - if I'm wrong and it's actually quite simple to fit 3 splines, can you show me how? Perhaps a good way to order my kind of input points would be by their distance from the start/end locations. I'm working on a function of the form
> function cPts = fitCurveTo3DPts(XYZ, stPt, endPt, n)
>
> I'll post here if/when it's done, but I'd appreciate any input you have.

Ok, so here's what I've come up with. I think it's reasonably concise and seems to be working well for my type of input. Currently I'm hard-coding 5 iterations based purely on my own testing. I imagine it could be generalised to exit based on some degree of stability in the resulting control points. I can certainly imagine the type of twisting or complicated input for which my function below would handle quite poorly. Any suggestions would be great.

Cheers,
Sven.


function cPts = fitCurveTo3DPts(XYZ, stPt, endPt, n)
MAX_ITERATIONS = 5;
completedIterations = 0;

% Initialise control pts linearly between start/end anchors
cPts = interp1(0:1, [stPt; endPt], linspace(0,1,n));

while completedIterations < MAX_ITERATIONS
    % Get the nearest-neighbour cntrl pt for each of the sample points
    sqDists = cellfun(@(x)sum(bsxfun(@minus, XYZ, x).^2,2), num2cell(cPts,2),'UniformOutput', false);
    [~, nnIdxs] = min(cat(2,sqDists{:}),[],2);
    % Keep the anchors, update inner cPts to the mean of their linked input pts
    for i = 2:n-1
        cPts(i,:) = mean(XYZ(nnIdxs==i,:),1);
    end
    % Handle any cPts that didn't have linked pts so their mean became NaN
    goodIdxs = find(~isnan(cPts(:,1)));
    badIdxs = find(isnan(cPts(:,1)));
    cPts(badIdxs,:) = interp1(goodIdxs, cPts(goodIdxs,:), badIdxs);
    % Re-spread the control points out linearly
    cPtCumSumDists = cumsum([0; sqrt(sum(diff(cPts,1).^2,2))]);
    cPts = interp1(cPtCumSumDists, cPts, linspace(0, cPtCumSumDists(end), n));
    completedIterations = completedIterations + 1;
end

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us