I may be wrong but it sounds like you are just trying to do a Total Least Squares fits to 3D data for different sets of data (lines)? This is a little tougher than doing regular 2D TLS, but you can do it in matlab pretty easily. An example using built-in matlab functions is at:
Alternatively, you could so it using linear algebra, by minimizing the ‘true’ distances between the points and the line. Representing your line as [x;y;z] = [a;b;c]*t + [x0;y0;z0], [a;b;c] is your normal vector along the line, and [x0;y0;z0] is a point on the line. You can take x0=mean(x), etc since a TLS fit will pass through the averages of the points. To find your normal vector, you can minimize the distances between the points and the points on the line perpendicular to them. For this setup, your system is a matrix representation of the cross products of your normal vector (n) and the vector from your [x0;y0;z0] point to your data points (lets call it r).
IE: ri = [xi,yi,zi] - [x0,y0,z0]
M = [n x r0;
n x ri
Alternatively, this minimization is equivalent to maximizing the projection of your vector r onto your normal vector, which is easier to setup.
IE: maximize the norm of A:
A = [x-mean(x), y-mean(y), z-mean(z)]
For either system, you end up taking the SVD of your system (A), and then your normal vector is the singular vector corresponding to either the maximum or minimum singular value, depending on which system you used.
note: I’m not familiar with the function used in the link I posted, but the ‘coeff’ variable they show in the example is the same as the singular vectors of the A matrix above.