Bounded high order polynomial regression
I am doing bone segmentation, whereas the result of this segmentation is the points located in a circular structure around this bone. However, as is customary with a qCT scan, there is quite a bit of noise (like flesh) at the points I have. So the general problem is how to remove this noise.
What I'm doing at the moment is converting it to polar coordinates of the shifted data, so that I get the data at the distance from the center of the bone to a point depending on the angle.
[THETA,RHO] = car2pol(N(1,:),N(2,:),center);
[THETA, id] = unique(THETA);
Once I have this data, I do a high order polynomial regression (30), which removes noise nicely and gives a nice smooth curve.
[p,~,mu] = polyfit(THETA,RHO,30); RHO = polyval(p,THETA,[],mu);
After that I convert it back to Cartesian coordinates:
[x,y] = pol2car(THETA,RHO,center);
My only problem here is that the start and end points are not necessarily the same point after the regression as they should be. So my question is, can you do some type of bounded polynomial regression where I can ensure that the y-value of the first point has the same value as the y-value of the other endpoint? Or is there any other way to do this?
source to share
I would recommend 4 point 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 considering the following points
>> plot(x, y, '.')
you can convert 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 do standard linear regression
>> b = regress(r, X);
Then you can get the predictions and convert back to 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)
source to share