Hi,
I'm sure this isn't the best solution, but if you enter f, Xa,Xb,Ya,Yb,Za, and Zb as symbols, and then enter your four constraints, and then your equation for D^2, and then use expand, you can pretty easily read off the terms and group them as you like in the quadratic form. I am sure there is a better way, but thought this might help you move along until someone has a better idea.
syms f Xa Xb Za Zb Ya Yb xa=-f*Xa/Za;
xb=-f*Xb/Zb;
ya=-f*Ya/Za;
yb=-f*Yb/Zb;
D2=((xa-xb)^2+(ya-yb)^2+(Za-Zb)^2);
expand(D2)
