Rank: 843 based on 163 downloads (last 30 days) and 4 files submitted
photo

Dan Couture

E-mail

Personal Profile:

https://github.com/MathYourLife/Matlab-Tools - Github Matlab tools
https://twitter.com/MathYourLife
www.MathYourLife.com

Professional Interests:

 

Watch this Author's files

 

Files Posted by Dan View all
Updated   File Tags Downloads
(last 30 days)
Comments Rating
26 Aug 2013 Screenshot Plane Fitting and Normal Calculation Given a set of x,y,z coordinates, find the best planar fit. Author: Dan Couture plot, regression, fit, modelling, plane, normal 99 24
  • 4.5
4.5 | 10 ratings
29 Aug 2012 Screenshot SYSLOG Utility with File Rotation A basic logging function with log rotation that records the message along with system the time Author: Dan Couture syslog, debug, logging, file 17 0
20 Aug 2012 Random String Utility Generate a random string. Specify length and type of characters. Author: Dan Couture random, string 19 2
13 Aug 2012 Screenshot Generate an Orthogonal Set of Unit Vectors Given a vector, a set of orthogonal unit vectors is calculated to use as rotated coordinate axes. Author: Dan Couture mathematics, plot, coordinate axes, orthogonal 28 0
Comments and Ratings by Dan View all
Updated File Comments Rating
26 Jun 2014 Plane Fitting and Normal Calculation Given a set of x,y,z coordinates, find the best planar fit. Author: Dan Couture

Good Morning John,

Agreed, additional context for error handling is beneficial. I'll add that in.

Thanks for the key word suggestions.

SVD vs linear regression:

Thanks for the high noise example illustrating the bias. The comparisons I run previously had relatively low noise ~ 0.1%. I used the following to compare results between SVD and linear regression:

function test_normals(N)
N = N/norm(N);
data = randn(1000,3);
data = data - (data*N)*N';
data = data + randn(size(data))/2;

n1 = fitNormal(data)';
n2 = svd_method(data)';

r_err = acos( abs(dot(n1, N)) / (norm(n1)*norm(n2)) ) / pi * 180
svd_err = acos( abs(dot(n2, N)) / (norm(n1)*norm(n2)) ) / pi * 180
svd_err / r_err

%plot3(data(:,1),data(:,2),data(:,3),'o')
%input("Press Enter")
end

for i = 1:10
N = [1 1 1]';
% or
N = rand(3,1);
% or
N = [0 0 1]';
test_normals(N)
end

For the worst case scenario for linear regressions N=[1 1 1]', SVD shows a significant improvement with ~ 1 degree error vs ~ 16 degree error for the regression. On the other hand, the best case scenario for regressions such as N=[0 0 1]' the regression techniques shows a modest improvement in accuracy over SVD ( ~1.1 degree error for regression vs ~1.3 degree error for svd )

Going back to the main reason I wrote this module. In the spectrum of brute force to elegant techniques SVD is clearly a more elegant and robust technique. For students and programmers that are early on in their mathematics or have been a while since they have broken out their TI's a more rudimentary technique using simple matrix operations may be easier to comprehend. I'm planning on incorporating many of the notes you suggested, but will stick with the regression method as I believe code sample should be available for multiple points along the spectrum of mathematics. I do like the idea of taking the opportunity to reference SVD and it's advantages over regression and will be adding notes to the description and comments to that point.

Dan

26 Jun 2014 Plane Fitting and Normal Calculation Given a set of x,y,z coordinates, find the best planar fit. Author: Dan Couture

Hi John,

Thanks for taking the time to put together such detailed notes. It seems like there were 5 main points.

Summary
------------------------------------------------------
1) H1 line
* good, but would prefer keywords

2) Error checking
* data array should be checked for Nx3
* show_graph should be checked for boolean type

3) In code comments and variable names
* approved

4) Check for determinant = 0
* Useless with floating point error
* should use cond() or rank()

5) Algorithm
* A/b is better than x = inv(A'*A)*A'*b
* SVD is the way it should be done

------------------------------------------------------
1) H1 line

Thank you. I'll look into adding some keywords. Any suggestions?

2) Error checking

I'll admit, allowing the code to throw an exception and expect it to be handled by the calling function is a personal preference, but I can see the merit in explicitly checking prior to use.

3) In code comments and variable names

Thanks again

4) Check for determinant = 0

I agree that this check is largely gratuitous as float values and noise will prevent the determinant from reaching zero. This check does prevent against errors in some edge cases though such as (if I can borrow from your sample)

N = [2 0 0]';
N = N/norm(N);
data = randn(10,3) * 1000;
data = data - (data*N)*N';

where there is no variation along the x axis and the determinant does evaluate to zero exactly.

I do like your suggestion of cond or rank vs det and will incorporate the change.

5) Algorithm

A cursory skim of your comment gave the impression that fitNormal was returning incorrect results. I compared results between the fitNormal and the SVD method you provided (below) and the returned vectors were equivalent.

mu = mean(data,1);
[U,S,V] = svd(bsxfun(@minus,data,mu));
normal_vector = V(:,3)

Do you have a sample data set where the results do not agree?

I believe what your main point is that the SVD method is a far more efficient method for calculating this this fit. To that, I wholeheartedly agree. The reason I stuck with the regression method, although less efficient, is that it only requires basic knowledge of matrix algebra (product, inverse, transpose, ...) that could be verified by hand by high school algebra students. Since I believe this module to be mainly used in exploration of data rather than in high performance situations I erred on the side of a less efficient method to allow a wider audience of coders that took the time to read the source to be able to verify the results themselves.

With that said, I am always in favor of exposing people to more sophisticated and rigorous mathematical techniques. I'll modify the comments in the header to clarify the circuitous approach used and refer them to the SVD method if greater efficiency is required.

Dan

25 Jun 2014 Plane Fitting and Normal Calculation Given a set of x,y,z coordinates, find the best planar fit. Author: Dan Couture

Hi Gaoyang,

That first `for i = 1:3` section generates the solutions for the least squares regression via matrix algebra for the three equations of x=f(y,z)+c, y=f(x,z)+c, and z=f(x,y)+c.

The standard way the equation is written is:
b =(X'X)-1 X'y

where X is a (n by m) matrix of input values and y is a (n by 1) vector.

Since the inverse of a matrix with a determinant of zero doesn't exist, checking for `det(X_m)==0` prevents an error being thrown at `coeff = (X_m)^-1 * X' * data(:,i);`

You are correct the regression generates a normal vector that describes the plane. Each of these vectors (for i=1,2,3) is later used to quantify the residuals from the data points and determine which is the best fit.

25 Aug 2013 Plane Fitting and Normal Calculation Given a set of x,y,z coordinates, find the best planar fit. Author: Dan Couture

No Problem. Thanks for the feedback. A new version was just submitted and should be accepted soon. If you need to you can also pull it from git https://github.com/MathYourLife/Matlab-Tools

22 Aug 2013 Plane Fitting and Normal Calculation Given a set of x,y,z coordinates, find the best planar fit. Author: Dan Couture

Thank for the feedback Amir. It's likely that your data doesn't fit in the window set by the code. It's hard coded now for [-1, 1] as you can see on this line.

set(get(L, 'Parent'),'DataAspectRatio',[1 1 1],'XLim',[-1 1],'YLim',[-1 1],'ZLim',[-1 1]);

I'll see if I can get that updated so that it jumps to wherever you data set is. I'll probably keep the range of all the axes equal though or it will make the plot appear as though the solution does not come out to be normal to the plane.

Comments and Ratings on Dan's Files View all
Updated File Comment by Comments Rating
27 Jun 2014 Plane Fitting and Normal Calculation Given a set of x,y,z coordinates, find the best planar fit. Author: Dan Couture Gaoyang

Still not sure why you change to negative, and assign the first element of the normal vectors.

% Construct and normalize the normal vector
coeff = (X_m)^-1 * X' * data(:,i);
c_neg = -coeff;
c_neg(i) = 1;
coeff(i) = 1;
n(:,i) = c_neg / norm(coeff);

27 Jun 2014 Plane Fitting and Normal Calculation Given a set of x,y,z coordinates, find the best planar fit. Author: Dan Couture Gaoyang

Hi Dan,

Thank you for your quick reply. Now I can understand your logic.

By the way, Niloy Mitra used the third principal component of point clouds as the surface normal. The logic is that the first two represent the biggest variation which is on the plane. and the third one is the surface normal.

27 Jun 2014 Plane Fitting and Normal Calculation Given a set of x,y,z coordinates, find the best planar fit. Author: Dan Couture Gaoyang

But I doubt that the robustness check for determinant is necessary. I used this before, due to the precision in Maltab, it will be very rare to have exact 'zero'. I also use some range to define 'zero', but it is data sensitive because that would be related the scaling (range) of data.

26 Jun 2014 Plane Fitting and Normal Calculation Given a set of x,y,z coordinates, find the best planar fit. Author: Dan Couture D'Errico, John

Yeah, with low noise, any method works reasonably well. With higher noise, the points on the high end that have error that brings them in have less leverage than those points that get moved out. This tends to introduce a bias.

Of course, it is your choice which method you choose. Definitely use backslash for any linear regressions instead of the normal equations and avoid det, and I will be happier.

You might offer both methods as an option, because whenever I make a choice that limits my code, I always find someone who would have preferred I made the alternative choice.

Anyway, I'm always willing to review submissions when they have been improved.

26 Jun 2014 Plane Fitting and Normal Calculation Given a set of x,y,z coordinates, find the best planar fit. Author: Dan Couture Couture, Dan

Good Morning John,

Agreed, additional context for error handling is beneficial. I'll add that in.

Thanks for the key word suggestions.

SVD vs linear regression:

Thanks for the high noise example illustrating the bias. The comparisons I run previously had relatively low noise ~ 0.1%. I used the following to compare results between SVD and linear regression:

function test_normals(N)
N = N/norm(N);
data = randn(1000,3);
data = data - (data*N)*N';
data = data + randn(size(data))/2;

n1 = fitNormal(data)';
n2 = svd_method(data)';

r_err = acos( abs(dot(n1, N)) / (norm(n1)*norm(n2)) ) / pi * 180
svd_err = acos( abs(dot(n2, N)) / (norm(n1)*norm(n2)) ) / pi * 180
svd_err / r_err

%plot3(data(:,1),data(:,2),data(:,3),'o')
%input("Press Enter")
end

for i = 1:10
N = [1 1 1]';
% or
N = rand(3,1);
% or
N = [0 0 1]';
test_normals(N)
end

For the worst case scenario for linear regressions N=[1 1 1]', SVD shows a significant improvement with ~ 1 degree error vs ~ 16 degree error for the regression. On the other hand, the best case scenario for regressions such as N=[0 0 1]' the regression techniques shows a modest improvement in accuracy over SVD ( ~1.1 degree error for regression vs ~1.3 degree error for svd )

Going back to the main reason I wrote this module. In the spectrum of brute force to elegant techniques SVD is clearly a more elegant and robust technique. For students and programmers that are early on in their mathematics or have been a while since they have broken out their TI's a more rudimentary technique using simple matrix operations may be easier to comprehend. I'm planning on incorporating many of the notes you suggested, but will stick with the regression method as I believe code sample should be available for multiple points along the spectrum of mathematics. I do like the idea of taking the opportunity to reference SVD and it's advantages over regression and will be adding notes to the description and comments to that point.

Dan

Contact us