I'm graphing some data I obtained in a physics experiment for school and I
have a set of points. I'm trying to find the best fit plane for the set of
points, I tried this algorithim from the book "3D Math Primer for Graphics
and Game Development" but I'm getting strange results and the plane is no
where near what it's supposed to look like.
All the points are already semi-planar so I'm assuming that the plane
shouldn't be orient nearly 90 degrees to the data. I'll post an example of
what going in in p.b.i if my explanation isn't very clear.

#local n = Number of points in array;

//Calculate Best Fit Normal
 #local Normal = <0, 0, 0>;
 #local P = Points[n-1];
 #local i = 0;
 #while(i < n)
  #local C = Points[i];
  #local Normal = Normal + <(P.z+C.z)*(P.y-C.y), (P.x+C.x)*(P.z-C.z),
(P.y+C.y)*(P.x-C.x)>;
  #local P = C;
  #local i = i + 1;
 #end
 #local Normal = vnormalize(Normal);

 //Calculate the best fit distance
 #local Distance = 0;
 #local i = 0;
 #while(i < n)
  #local Distance = Distance + vdot(Points[i], Normal);
  #local i = i + 1;
 #end
 #local Distance = Distance / n;

plane {
    Normal, Distance
}

any help with this issue would be greatly appreciated

--
Kevin
http://www.geocities.com/qsquared_1999/
#macro _(r)#if(r<12)#local i=asc(substr("oqshilacefg"
,r,1))-97;disc{<mod(i,7)-3,div(i,7)-1,6>,z,.4
pigment{rgb 10}}_(r+1)#end#end _(1)//KL

Post a reply to this message

Wasn't it Kevin Loney who wrote:
>I'm graphing some data I obtained in a physics experiment for school and I
>have a set of points. I'm trying to find the best fit plane for the set of
>points, I tried this algorithim from the book "3D Math Primer for Graphics
>and Game Development" but I'm getting strange results and the plane is no
>where near what it's supposed to look like.
>All the points are already semi-planar so I'm assuming that the plane
>shouldn't be orient nearly 90 degrees to the data. I'll post an example of
>what going in in p.b.i if my explanation isn't very clear.

The code you posted works perfectly for me. My guess would be that you
have a silly mistake somewhere in the code that calls the bit that you
posted. (Wild guess: perhaps it's a left/right handed co-ordinate system
problem)

Here's an example of a complete scene using the plane-fitting code that
you posted. The 20 points are displaced from the plane "y=x+0.5z-1" by
random amounts. I reckon that the plane produced by the best fit
algorithm is a pretty good fit to the data.


global_settings {assumed_gamma 1.0}

camera {location  <0,0,-10> look_at   <0.5,-0.3,0>  angle 12}

light_source { <-30, 100, -30> color rgb 1}

#declare Points = array[20];
  
#declare R=seed(2);
#declare i=0;
#while (i<20)
  #declare X=rand(R);
  #declare Z=rand(R);
  #declare Y=(X + 0.5*Z - 1 -rand(R)*0.2);
  #declare Points[i] = <X,Y,Z>;
  sphere {Points[i],0.04 pigment {rgb x}}
  #declare i=i+1;
#end

#local n = 20;

//Calculate Best Fit Normal
 #local Normal = <0, 0, 0>;
 #local P = Points[n-1];
 #local i = 0;
 #while(i < n)
  #local C = Points[i];
  #local Normal = Normal + <(P.z+C.z)*(P.y-C.y), (P.x+C.x)*(P.z-C.z),
(P.y+C.y)*(P.x-C.x)>;
  #local P = C;
  #local i = i + 1;
 #end
 #local Normal = vnormalize(Normal);

 //Calculate the best fit distance
 #local Distance = 0;
 #local i = 0;
 #while(i < n)
  #local Distance = Distance + vdot(Points[i], Normal);
  #local i = i + 1;
 #end
 #local Distance = Distance / n;

plane {Normal, Distance pigment {rgbt <1,1,1,0.7>}}

Post a reply to this message

Kevin Loney <klo### [at] pt2mcom> wrote:
>   #local Normal = Normal + <(P.z+C.z)*(P.y-C.y), (P.x+C.x)*(P.z-C.z),
> (P.y+C.y)*(P.x-C.x)>;

  I think that there might be a problem with the '+' operations there
(don't they need to be '-'?). To avoid any typing mistakes and to
optimize/simplify the code, you can write it like this:

#local Normal = Normal + vcross(P,C);

  However, what I wonder is if this can work. Some of the normals will
point to one side of the "plane" and others will point to the other.
Their sum can point virtually anywhere, depending on the arrangement of
the points.
  If the points were sorted in clockwise (or the anti-clockwise) order,
then it will probably ok.

-- 
#macro N(D)#if(D>99)cylinder{M()#local D=div(D,104);M().5,2pigment{rgb M()}}
N(D)#end#end#macro M()<mod(D,13)-6mod(div(D,13)8)-3,10>#end blob{
N(11117333955)N(4254934330)N(3900569407)N(7382340)N(3358)N(970)}//  - Warp -

Post a reply to this message

Warp <war### [at] tagpovrayorg> wrote in news:3e69f02c@news.povray.org:

> Kevin Loney <klo### [at] pt2mcom> wrote:
>>   #local Normal = Normal + <(P.z+C.z)*(P.y-C.y), (P.x+C.x)*(P.z-C.z),
>> (P.y+C.y)*(P.x-C.x)>;
>   I think that there might be a problem with the '+' operations there
> (don't they need to be '-'?).

No. It works as it is.

Note that this:
<
  (P.z - C.z)*(P.y - C.y),
  (P.x - C.x)*(P.z - C.z),
  (P.y - C.y)*(P.x - C.x)
>
is NOT equal to vcross(P, C)


> To avoid any typing mistakes and to
> optimize/simplify the code, you can write it like this:
> 
> #local Normal = Normal + vcross(P,C);
> 
>   However, what I wonder is if this can work. 
...

It won't.


But this should:

#declare Normal = Normal + vcross(P - C, P + C);



Tor Olav

Post a reply to this message

Tor Olav Kristensen <tor_olav_kCURLYAhotmail.com> wrote in 
news:Xns### [at] 204213191226:

> Warp <war### [at] tagpovrayorg> wrote in news:3e69f02c@news.povray.org:
...
>> To avoid any typing mistakes and to
>> optimize/simplify the code, you can write it like this:
>> 
>> #local Normal = Normal + vcross(P,C);
>> 
>>   However, what I wonder is if this can work. 
> ...
> 
> It won't.

Sorry Warp. I think your suggestion will work as well.


Tor Olav

Post a reply to this message

your right it does work, but for some weird and obscure reason not with any
of the sets of data I'm using :-P So I'm doing it the long way, find the
plane passing through every possible combination of three points and average
them, it's not the most accurate but neither is the data I'm working with.

Thanks for the help

--
Kevin
http://www.geocities.com/qsquared_1999/
#macro _(r)#if(r<12)#local i=asc(substr("oqshilacefg"
,r,1))-97;disc{<mod(i,7)-3,div(i,7)-1,6>,z,.4
pigment{rgb 10}}_(r+1)#end#end _(1)//KL

Post a reply to this message