Translation of a 3D plane
Afficher commentaires plus anciens
I have a set of 3D points out of which I have fitted a plane using Principal Component Analysis. Now, I want to generate/create a plane parallel to this one at a particular distance or passing through a point and find how these individual points project on that new plane.
Any help or input is appreciated.
The following code is an example.
randn('state',0);
X = mvnrnd([0 0 0], [1 .2 .7; .2 1 0; .7 0 1],50);
plot3(X(:,1),X(:,2),X(:,3),'bo');
grid on;
maxlim = max(abs(X(:)))*1.1;
axis([-maxlim maxlim -maxlim maxlim -maxlim maxlim]);
axis square
view(-23.5,5);
[coeff,score,roots] = princomp(X);
basis = coeff(:,1:2)
normal = coeff(:,3)
n,p] = size(X);
meanX = mean(X,1);
Xfit = repmat(meanX,n,1) + score(:,1:2)*coeff(:,1:2)';
residuals = X - Xfit;
error = abs((X - repmat(meanX,n,1))*normal);
sse = sum(error.^2)
[xgrid,ygrid] = meshgrid(linspace(min(X(:,1)),max(X(:,1)),5), ...
linspace(min(X(:,2)),max(X(:,2)),5));
zgrid = (1/normal(3)) .* (meanX*normal - (xgrid.*normal(1) + ygrid.*normal(2)));
h = mesh(xgrid,ygrid,zgrid,'EdgeColor',[0 0 0],'FaceAlpha',0);
hold on
above = (X-repmat(meanX,n,1))*normal > 0;
below = ~above;
nabove = sum(above);
X1 = [X(above,1) Xfit(above,1) nan*ones(nabove,1)];
X2 = [X(above,2) Xfit(above,2) nan*ones(nabove,1)];
X3 = [X(above,3) Xfit(above,3) nan*ones(nabove,1)];
plot3(X1',X2',X3','-', X(above,1),X(above,2),X(above,3),'o', 'Color',[0 .7 0]);
nbelow = sum(below);
X1 = [X(below,1) Xfit(below,1) nan*ones(nbelow,1)];
X2 = [X(below,2) Xfit(below,2) nan*ones(nbelow,1)];
X3 = [X(below,3) Xfit(below,3) nan*ones(nbelow,1)];
plot3(X1',X2',X3','-', X(below,1),X(below,2),X(below,3),'o', 'Color',[1 0 0]);
hold off
maxlim = max(abs(X(:)))*1.1;
axis([-maxlim maxlim -maxlim maxlim -maxlim maxlim]);
axis square
view(-23.5,5);
4 commentaires
David Young
le 5 Avr 2011
It's not possible to read your code. Please edit your question, select the code, and click the "{}Code" button, so that it is formatted. Thanks.
Phanikumar Bhamidipati
le 5 Avr 2011
Andrew Newell
le 5 Avr 2011
Is this code not doing what you want?
Phanikumar Bhamidipati
le 7 Avr 2011
Réponses (1)
Tim Zaman
le 6 Avr 2011
0 votes
Catégories
En savoir plus sur Kernel Distribution dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!