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!

"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) ),

"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?

"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..

"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

"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

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.

"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

"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];

"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