Constrained high order polynomial regression

Question

I am doing some bone segmentation whereas the result of this segmentation is points placed in a circular pattern around this bone. However as it is taken using a qCT scan, there is quite a lot noise (from e.g. flesh) on the points that i have. So the overall problem is how to remove this noise.

What I am doing at the moment is to transform it into polar coordinates on the shifted data, such that I get data on the distance from the center of the bone to the point depending on the angle.

[THETA,RHO] = car2pol(N(1,:),N(2,:),center);
[THETA, id] = unique(THETA);

Once I have those data I make a high order (30) polynomial regression, which removes the noise nicely and yields a nice smooth curve.

[p,~,mu] = polyfit(THETA,RHO,30);
RHO = polyval(p,THETA,[],mu); 

Whereafter I transform it back into cartesian coordinates:

[x,y] = pol2car(THETA,RHO,center);

My only problem here is the start and end point are not necessarily the same point after the regression, which they should be. So my question is, can you make some type of constrained polynomial regression, where I can enforce that the y-value of the first point have have the same value as the y-value of the other end point? Or are there any other way to do this?

Answer 1

I would recommend a fourier regression, rather than polynomial regression, i.e.

rho = a0 + a1 * cos(theta) + a2 * cos(2*theta) + a3 * cos(3*theta) + ...
           b1 * sin(theta) + b2 * sin(2*theta) + b3 * sin(3*theta) + ...

for example, given the following points

>> plot(x, y, '.')

you can transform to polar coordinates

>> [theta, r] = cart2pol(x, y);

and create an array of sines and cosines

>> X = ones(size(theta, 1), 1);
>> for n = 1:N
       X(:, end+1) = cos(n * theta);
       X(:, end+1) = sin(n * theta);
   end

and run a standard linear regression

>> b = regress(r, X);

You can then get the forecasts and transform back into cartesian coordinates to plot it

>> rhat = X * b;
>> [xhat,yhat] = pol2cart(theta, rhat);

which gives

>> plot(x, y, '.');
>> hold on;
>> plot(xhat, yhat, 'r', 'LineWidth', 2)