question about code for fitting plane to points

[cosinekitty]

Hi, I was very pleased to find your code for fitting a plane to a collection of possibly noisy points. I'm adapting this for some code I'm writing to analyze planetary orbits. There is something about the code at

https://www.ilikebigbits.com/2017_09_25_plane_from_points_2.html

that I don't understand. This part seems suspicious to me...

    let mut weighted_dir = Vec3{x: 0.0, y: 0.0, z: 0.0};

    {
        let det_x = yy*zz - yz*yz;
        let axis_dir = Vec3{
            x: det_x,
            y: xz*yz - xy*zz,
            z: xy*yz - xz*yy,
        };
        let mut weight = det_x * det_x;
        if weighted_dir.dot(&axis_dir) < 0.0 { weight = -weight; }
        weighted_dir += &axis_dir * weight;
    }

When you take the dot product of weighted_dir with axis_dir, the result is always going to be zero, so weight will always be left non-negative. I know what you're probably thinking... it sounds like I'm going to complain about efficiency.

Actually, the part that concerns me is that it looks like the code will produce different results based on the order you process the axes. Let's call the block of code after weighted_dir is initialized to (0,0,0) the "det_x block". There is also a det_y block and a det_z block below it. If I were to swap the det_x and det_y blocks, for example, now the det_x block will sometimes execute { weight = -weight; } (because now weighted_dir is no longer a zero vector) when before it would never do that.

It seems suspicious that the results depend on the order you add terms to the weighted_dir vector. Intuitively, the axes should all be treated equally and independently. Am I missing something?

Thanks for publishing this. Your writing is very clear and helps me understand the matrix stuff a lot better!